AD-A174  612  THE  FORHARD  AND  ADJOINT  TRANSPOST  OF  BULLETS  IN  HONTfe  1/1 
CARLO  AIR-OEFENSE-E  <U>  ARMY  BALLISTIC  RESEARCH  LAB 
ABERDEEN  PROVING  GROUND  MD  U  B  BEVERLY  JUL  86 
UNCLASSIFIED  BRL-TS-2758  F/G  19/6  NL 


AD- A 174  612 


us  ARMY 


MATERIEL 

COMMAND 

TECHNICAL  REPORT  BRL-TR-2750 


THE  FORWARD  AND  ADJOINT  TRANSPORT  OF 
BULLETS  IN  MONTE  CARLO  AIR-DEFENSE- 

END-GAME  BALLISTIC 
VULNERABILITY  PROBLEMS 


William  B.  Beverly 

Q_ 

8  July  1986 

LxJ 

fcd 
»  1 1 


APPROVED  FDR  PUMJC  RELEASE:  WSIRWHIOM  UNUMIIED. 


US  ARMY  BALLISTIC  RESEARCH  LABORATORY 

ABERDEEN  PROVING  GROUND,  MARYLAND 


DTIC 

ELECTE M 

la,  nco  n  irtoc  J 


%.  DEC  2 


1986 


Destroy  this  report  when  it  is  no  longer  needed. 
Do  not  return  it  to  the  originator. 


Additional  copies  of  this  report  may  be  obtained 
from  the  National  Technical  Information  Service, 

U.  S.  Department  of  Commerce,  Springfield,  Virginia 

22161. 


The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position,  unless  so  designated  by  other 
authorized  documents. 


The  use  of  trade  names  or  manufacturers'  names  in  this  report 
does  not  constitute  indorsement  of  any  commercial  product. 


UNCLASSIFIED _ 

SECURITY  CLASS!  Fl  CATION  OF  THIS  PAGE  (Whan  Dmtm  Entered) 

1  _  REPORT  DOCUMENTATION  PAGE 


_ REPORT  DOCUMENTATION  PAGE _  bm^™™tS^foi«< 

|T  REPORT  NUMBER  |2.  GOVT  ACCESSION  NO.  ■*  RECIPIENT'S  CATALOG  NUMBER 

Technical  Report  BRL-TR-2750  Ab- A  1*1^  &  f  '  ' _ _ 

4.  TITLE  fand  Subtitle)  i-  TYPE  OF  REPORT  4  PERIOD  COVER 

The  Forward  and  Adjoint  Transport  of  Bullets  in 
Monte  Carlo  Air-Defense-End-Game  Ballistic 

Vulnerability  Problems  6.  performing  org.  report  number 


:  s.  type  of  report  &  period  covered 


I  6.  PERFORMING  ORC.  REPORT  NUMBER 


7.  authorc*; 

William  B.  Beverly 


8.  CONTRACT  OR  GRANT  NUMBERr«> 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS  <0-  PROGRAM  ELEMENT.  PROJECT.  TASK 

US  Army  Ballistic  Research  Laboratory  area  «  work  unit  numbers 

ATTN:  SLCBR-SE 

Aberdeen  Proving  Ground,  MD  21005-5066 

II.  CONTROLLING  OFFICE  NAME  ANO  ADORESS  12.  REPORT  DATE 

US  Army  Ballistic  Research  Laboratory  Jul  2.986 

ATTN:  SLCBR-DD— T  u.  number  of  pages 

Aberdeen  Proving  Ground,  MD  21005-5066  72 

14.  MONITORING  AGENCY  NAME  8  ADDRESS^//  dllletent  from  Controlling  Olllce)  IS  SECURITY  CLASS.  ( ot  Ihle  report) 


Unclassified 


15a.  DECL  ASSIFl  CATION /DOWNGRADING 
SCHEDULE 


16.  DISTRIBUTION  STATEMENT  (o I  Ihl a  Report) 

Approved  for  public  release.  Distribution  unlimited. 


[  17.  DISTRIBUTION  STATEMENT  (ol  the  abstract  entered  In  Block  20,  II  dlllerenl  Iron  Report) 


18.  SUPPLEMENTARY  NOTES 


19.  KEY  WORDS  (Continue  on  r«r«r««  old*  It  n«c«a«ary  and  Identity  by  block  numbor) 

Monte  Carlo  Methodology 
Air  Defense  Calculations 
Expected  Value 
Kill  Probabilities 

ZO.  ABSTRACT  fC aertBeue  an  rerereo  «S  II  neceeeery  mrd  Identity  by  block  number)  (  i  rik  ) 

Computerized  Monte  Carolo  methodologies  are  routinely  used  to  estimate  the 
effectiveness  of  weapon  systems.  In  particular,  such  codes  are  used  at  this 
installation  to  simulate  the  performance  of  air  defense  gun  systems.  In  gen¬ 
eral,  these  methodologies  are  constructed  to  conduct  sample  firing  events  by 
picking  values  for  the  system  uncertainties  from  their  natural  probability  den¬ 
sity  functions.  Therefore,  such  methodologies  can  be  easily  (Continued  on 
next  page) 


1473  EDITtOM  OF  *  NOV  69  IS  OBSOLETE 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PA'.E  f When  Dele  Fnrcr.tJ) 


UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  of  THIS  PAGE{TOi«n  Dmlm  Enltrmd) 


constructed  with  a  logic  which  compares  one-to-one  with  real-life  firing  events 
without  ever  constructing  an  explicit  formulation  of  the  associated  expected- 
value-of-kill  integrals.  However,  this  ease  of  construction  is  not  an  unmiti¬ 
gated  advantage  since  an  analytical  formulation  may  permit  the  development  of 
more  efficient  procedures  for  calculating  conventional  kill  probabilities  as 
well  as  providing  the  basis  for  conducting  more  elaborate  studies  such  as  error 
and  parameter-sensitivity  studies.  For  example,  in  a  "natural"  Monte  Carlo 
simulation,  the  fractional  uncertainty  in  a  kill  estimate  increases  as  an  in¬ 
creasing  fraction  of  the  sample  bullets  miss  the  target.  The  efficiency  of  such 
a  calculation  may  be  improved  if  procedures  could  be  developed  where  most  of 
the  sample  bullets  hit  the  target. 

Therefore,  our  first  objective  in  this  report  is  to  develop  a  mathematical 
formulation  whose  modeling  of  weapon-system  uncertainties  is  similar  to  the  cur^ 
rent  methodologies  being  used.  To  accomplish  this,  we  formulate  the  expected- 
value-of-kill  integral  where  the  variables  of  integration  are  those  used  in  the 
derivation  of  a  forward,  or  "natural",  Monte  Carlo  evaluation.  Then,  in  an 
attempt  to  develop  alternate  formulation  which  may  permit  the  more  efficient 
calculation  of  kill  probabilities,  we  transform  the  integrand  to  an  adjoint 
form  such  that,  in  its  Monte  Carlo  evaluation,  the  majority  of  the  sample  bullet 
trajectories  will  intersect  the  target.  We  then  outline  general  procedures  for 
conducting  Monte  Carlo  evaluations  of  both  the  original  and  the  transformed 
expected-value-of-kill  integrals.  Finally,  we  devise  an  example,  two- 
dimensional  problem,  formulate  the  expected-value-of-kill  integral  to  both  the 
"natural"  and  the  adjoint  mode,  and  compare  the  efficiencies  of  the  Monte  Carlo 
evaluation  of  each  mode  for  a  series  of  different  scenarios. 


TABLE  OF  CONTENTS 


Page  No. 


LIST  OF  FIGURES . 5 

I.  INTRODUCTION . 7 

II.  THE  CONSTRUCTION  OF  EXPECTED-VALUE-OF-KILL 

INTEGRALS . 16 

III.  A  GENERAL  MONTE  CARLO  EVALUATION  OF  THE 

NATURAL  AND  THE  TRANSFORMED  EXPECTED- 
VALUE-OF-KILL  INTEGRALS . 27 

IV.  THE  SOLUTION  OF  A  PEDAGOGICAL  EXAMPLE  AIR-DEFENSE-END¬ 

GAME  PROBLEM iUSING  FORWARD  AND  ADJOINT  SAMPLING 
PROCEDURES . 3  6 


A.  The  Problem  Set  Description . 36 

B.  The  Forward  Monte  Carlo  Estimation  Procedures . 42 

C.  The  Adjoint  Monte  Carlo  Estimation  Procedures . 46 

D.  The  Problem  Results . 54 

V.  CONCLUSIONS  AND  PROGNOSIS . 57 

REFERENCES . 60 

APPENDIX  . ttt, . 63 

DISTRIBUTION  LIST . 71 


-3- 


LIST  OF  FIGURES 


Figure  No.  Page  No. 

1.  The  Line-of-Sight  (LOS)  Range  Vector  for  an  Example 
Air-Defense  End-Game  Bullet  Trajectory  (No  Wind 

Effects) . 9 

2.  The  Vectors  (Gun-Coordinate  System)  Associated  with 

an  Air-Defense  End-Game  Problem . 14 

3.  An  Illustration  of  a  Bullet-Trajectory  Impact  with  an 

Example  Surface  S( . 17 

4.  The  Target  Spherical  Coordinate  System  Associated 

with  the  Change  of  Integration  Variables . 23 

5.  The  Orientation  of  the  Gun-Coordinate  System  and 
the  Target-Coordinate  System  When  the  Bullet 

Intersects  S(\rb  =  (rt)(] . 24 

6.  The  Forward,  or  “Natural”,  Monte  Carlo  Estimation  of 
the  Expected-Value-of-Kill  of  a  Critical  Component  in 

a  Target  Aircraft . 28 

7.  The  Adjoint  Monte  Carlo  Estimation  of  the  Expected- 
Value-of-Kill  of  a  Critical  Component  in  a  Target 

Aircraft . 31 

8.  A  Two-Dimensional  Air-Defense  End-Game  Example 

Problem  (Not  To  Scale) . 37 

9.  The  Gun-Coordinate  System  (GCS)  for  the  Example 

Problem  (Not  To  Scale) . 38 

10.  The  Target-Coordinate  System  (TCS)  for  the  Example 

Problem  (Not  To  Scale) . 39 

11,  The  Probability  Density  Functions  Associated  with  the 

Example  Problem . 40 


-5- 


LIST  OF  FIGURES  (Continued) 


Figure  No.  Page  No. 

12.  The  Forward  Monte  Carlo  Estimation  of  the 

Expected-Value-of-Hit  for  the  Example  Problem . 43 

13.  The  Time-of-Flight  Analytic  Right  Triangle  Used 

Used  To  Calculate  Aiming  Angles  (Not  To  Scale) . 45 

14.  The  Time-of-Flight  Analytic  Right  Triangle  Which 
Relates  Target  Hits  in  the  GCS  and  the  TCS  (Not 

To  Scale) . 47 

15.  The  Adjoint  Monte  Carlo  Estimation  of  the  Expected- 

Value-of-Hit  for  the  Example  Problem . 49 

16.  An  Illustration  of  the  Limits  Associated  with  the 

Target  Polar  Angle  <f>  (Not  To  Scale) . 52 

17.  A  Graphic  Comparison  of  the  Calculational  Efficiencies 

of  the  Forward  and  Adjoint  Monte  Carlo  Procedures . 56 


-6- 


I.  INTRODUCTION 


Ongoing  work  at  the  US  Army  Ballistic  Research  Laboratory  (BRL)  is 
aimed  at  predicting  the  effectiveness  of  air-defense  gun  systems.  Computer¬ 
ized  methodologies  are  used  to  solve  representative  end  game  problems. 
One  such  methodology,  illustrated  by  the  MGEM^  Monte  Carlo  computer 
program, r  is  currently  used  at  BRL  to  calculate  the  expected  value  of  kill 
for  aircraft  exposed  to  point  detonation r(PD)  bullets  fired  from  gun  systems 
located  on  the  ground.  However,  the  procedures  used  in  MGEM  are  not 
adequate  for  predicting  the  effectiveness  of  proximity  fuzed  (PX)  bullets 
and  need  to  be  extended  to  obtain  a  methodology  with  such  a  capability. 
Additionally,  the  need  for  even  more  complex  methodologies  is  anticipated 
in  the  future  when  the  feasibility  of  using  "smart”  PX  bullets,  that  is  bul¬ 
lets  whose  trajectories  can  be  changed  during  flight  in  order  to  take  advan¬ 
tage  of  radar  fixes  of  target  location  at  times  subsequent  to  the  firing  of  the 
bullet,  is  investigated.  The  objectives  of  this  study  are  to  formalize  the 
MGEM  methodology,  and  to  develop  alternative  Monte  Carlo  methods  of 
solution  which  may  be  advantageously  used  in  the  future  in  the  construc¬ 
tion  of  PX  and  "smart-bullet^methodologies.  “  "  '■ 

The  expected  effectiveness  of  a  gun  system  depends  upon  the  accumula¬ 
tion  of  the  uncertainties  associated  with  target  detection  and  prediction  of 
its  future  locations,  the  prediction  of  the  bullet  flight  toward  a  target  inter¬ 
ception,  and  the  incapacitation  of  the  aircraft  by  a  hit.  For  the  purpose  of 
our  discussion,  we  have  lumped  these  uncertainties  into  nine  general 
categories  which  are  tabulated  below.  In  this  tabulation  and  in  the  follow¬ 
ing  discussion,  we  will  use  the  standard  conventions  of  probability  as 
described  by  C.  Eisenhart  and  M.  Zelen2  whereby  a  random  variable  is 
identified  by  using  an  upper  case  letter  and  its  values  are  identified  by 
using  the  associated  lower  case  letter. 

In  no  particular  order,  the  uncertainties  and  the  random  variables  used 
to  define  these  uncertainties  along  with  other  necessary  variables  and  nota- 
tional  conventions  are  tabulated  below: 


,R.A.  Scheder,  “Modern  Fire  Control  Analysis  Using  the  Modern  Gun  Effectiveness 
Model  (MGEM),”  U.S.  Army  Materiel  Systems  Analysis  Activity  Technical  Report  No. 
180,  October  1976. 

*C.  Eisenhart,  and  M.  Zelen,  “Elements  of  Probability,”  Handbook  of  Physics.  E.U. 
Condon,  and  H.  Odishaw,  Editors,  McGraw-Hill  Book  Cor  >any,  Inc.,  New  York,  1958. 
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1.  Variation  between  Rounds.  One  result  of  unavoidable  differences 
between  individual  rounds  of  ammunition  is  a  spread  in  the  muzzle  veloci¬ 
ties  Vm  of  the  associated  bullets. 

2.  Bore-Sight  Error  in  the  Gun.  This  error  causes  an  uncertainty  in 
the  direction  (elevation  angle  ©0  and  azimuth  angle  ♦(,)  of  the  bullet  as  it 
leaves  the  gun.  The  subscript  0  is  used  here  and  elsewhere  to  identify  vari¬ 
ables  associated  with  the  birth  of  a  bullet,  taken  to  be  the  muzzle  of  the 
gun.  Associated  with  these  firing  angles,  we  will  also  identify  here  the  line- 
of-sight  (LOS)  range  /  which  is  taken  to  be  the  distance  ai  some  flight  time 
tj  that  a  bullet  would  travel  in  still  air  with  no  gravity  (Figure  l)  The 
values  of  (0O >0)  are  assembled  with  l(tj)  to  obtain  the  LOS  range  vector 
[0o>V’o>4t/)l  which  is  used  to  calculate  the  location  of  a  bullet  during  its  flyout. 

The  reader  in  referring  to  Reference  1  should  note  that  our  “/”  is  the 
same  as  his  “r”.  We  have  reserved  the  bold  face  letter  “r"  with  an 
appropriate  subscript  to  identify  position  vectors,  while  as  noted  at  the 
time  of  use,  the  non-bold  face  letter  “r”  is  the  range  of  the  bullet  in  the 
spherical  coordinate  system  of  the  gun. 

Depending  on  analytical  convenience,  a  bullet  can  be  located  in  either  a 
Gun  Coordinate  System  (GCS)  or  a  Target  Coordinate  System  (TCS).  The 
GCS  is  located  with  its  origin  at  the  center  of  the  gun  muzzle  which  is  also 
taken  to  be  the  birth  site  of  bullets.  Depending  on  the  analytical  conveni¬ 
ence  of  the  moment,  either  spherical  coordinates  (0,ip,r)  or  Cartesian  coordi¬ 
nates  (x,y,z)  will  be  used.  The  polar  angle  6  is  measured  from  the  positive  z- 
axis  and  the  azimuthal  angle  tp  is  measured  from  the  positive  x-axis.  In 
compact  notation,  a  bullet  is  located  in  the  GCS  by  the  position  vector  rb 
and  the  reference  point  in  the  target  aircraft  is  located  by  the  position  vec¬ 
tor  *v  The  TCS  will  be  defined  when  needed  in  Section  II  of  this  report. 

3.  Dispersion  between  Shots  in  a  Burst.  The  servo  controlled  return  of 
the  gun  to  a  ready-to-fire  state  between  shots  in  a  burst  is  not  precise. 
Consequently,  an  additional  uncertainty  in  the  initial  direction  of  motion  of 
a  bullet  is  added  to  Item  2  where  the  uncertainty  is  expected  to  increase 
with  successive  firings  in  a  burst. 

4.  Variation  in  Air  Resistance  to  the  Flight  of  a  Bullet.  Different 
orientations  of  a  bullet  with  respect  to  its  trajectory  will  cause  the  resis¬ 
tance  by  air  to  vary.  This  effect  is  quantified  as  an  uncertainty  in  the  drag 
coefficient  Cd.  In  both  MGEM  and  the  present  stojdy,  the  uncertainty  is 
taken  to  be  Gaussian  where  Cd  is  assumed  to  have  the  same  value  over  the 
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Associated  Bullet 
Trajectory 


Figure  1.  The  Line-of-Sight  (LOS)  Range  Vector  for  an  Example  Air- 
Defense  End-Game  Bullet  Trajectory  (No  Wind  Effects) 


Mv.v: 


total  trajectory  of  a  bullet. 

5.  Variations  in  Wind.  The  effect  of  wind  on  a  trajectory  may  differ 
from  the  wind  corrections  used  to  calculate  aiming  angles.  This  efTeet,  ( W , 
and  \Vy),  is  assumed  to  be  appreciable  only  in  the  horizontal  plane.  In  both 
MGEM  and  the  present  study,  this  uncertainty  is  also  taken  to  be  Gaussian 
where  W.  and  Wu  are  taken  to  have  the  same  sets  of  values  over  the  total 
trajectory  of  a  bullet. 

6.  Time  of  Firing.  We  allow  in  our  formulation  for  an  uncertainty  in 
the  time  of  firing  T0  although  this  random  variable  is  assumed  to  be  pre¬ 
cisely  known  in  the  MGEM  methodology. 

7.  Target  Maneuver.  The  path  of  a  target  aircraft  from  the  time  of 
firing  of  a  bullet  to  its  anticipated  interception  by  the  bullet  can  not  be 
predicted  by  the  gun  system  with  absolute  certainty.  In  common  with 
MGEM,  we  will  assume  that  some  method  of  optimal  estimation  such  as 
the  Kalman  filter3  is  used  to  predict  the  most  probable  point  of  intercep¬ 
tion  derived  from  earlier  radar  fixes  of  target  location  and  velocity.  In 
Monte  Carlo  calculations,  a  target  maneuver  is  selected  from  a  collection  of 
standard  maneuver  paths  where  each  path  is  composed  of  the  aircraft  posi¬ 
tion,  velocity,  acceleration,  and  orientation  at  regular  time  intervals,  usu¬ 
ally  0.05  seconds. 

8.  Errors  in  the  Radar  Fixing  of  a  Target  Location.  The  true  location 
and  velocity  of  a  target  aircraft  may  differ  from  that  determined  by  a  fire- 
control  radar.  This  will  lead  to  an  additional  uncertainty  in  the  predicted 
future  locations  of  the  target  with  a  resultant  uncertainty  in  the  calculated 
aiming  angles  for  the  gun.  We  identify  the  estimated  aiming  angles  which 
are  derived  from  an  optimal  estimation  of  future  target  locations  based  on 
a  particular  set  of  radar  fixes  as  (0O,*O).  The  values  for  this  pair  of  random 
variables  are  taken  to  be  the  mean  values  of  the  firing  angles  described  ear¬ 
lier  in  Item  2.  In  the  present  study,  we  assume  that  all  of  this  uncertainty 
is  acquired  during  the  detection  and  analysis  of  the  radar  signals,  and  that 
the  relative  location  of  the  radar  in  the  GCS  is  precisely  known. 


3A.  Gelb,  J.F.  Kasper,  Jr.,  R.A.  Nash,  Jr.,  C.F.  Price,  and  A. A.  Sutherland,  Applied 
Optimal  Estimation.  A.  Gelb,  Editor,  The  M.I.T.  Press,  Massachusetts  Institute  of 
Technology,  Cambridge,  Massachusetts,  and  London,  England,  1974. 
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9.  The  Kill  Probability  of  a  Critical  Component  in  the  Target  Aircraft 
by  a  Bullet.  In  general,  the  causal  prediction  of  some  specified  loss  of  capa¬ 
bility  by  a  component  due  to  a  bullet  is  a  huge  task  which  would  require  an 
impractical,  if  not  impossible,  detailed  calculation  of  the  bullet-target 
interaction.  The  description  of  all  aspects  of  this  phenomena  would  require 
a  large  number  of  random  variables;  for  example,  the  generation  of  metal 
fragments  from  a  detonating  PX  round  and  the  penetration  of  these  metal 
fragments  into  the  target  aircraft  may  require  a  number  of  random  vari¬ 
ables  approaching  Avagodro’s  number  since  the  crystalline  structure  of  the 
bullet  and  the  target  are  involved.  Instead,  attempts  have  been  made  in 
ballistic-vulnerability  methodologies  to  characterize  the  behavior  of  a  large 
aggregate  of  similar  rounds  of  a  type  by  a  small,  manageable  number  of 
random  variables  and  hope  that  the  true  expeeted-value-of-capability  loss  is 
adequately  approximated  by  a  function  of  these  used  random  variables. 
For  example,  according  to  Kruse  and  Brizzolara4,  the  P K\H  (probability  of  a 
kill  given  a  hit)  function  for  a  component  is  generally  given  in  terms  of  the 
mass  and  velocity  of  a  metal  fragment  as  it  impacts  the  outer  surface  of  a 
critical  component.  An  attempt  to  improve  the  mass-velocity  characteriza¬ 
tion  of  metal  fragments  is  also  mentioned  in  this  reference. 

Item  9  completes  the  tabulation  of  the  lumped  uncertainties  associated 
with  air-defense  end-game  problems.  In  general,  the  probability  density  of 
random-variable  values  associated  with  an  uncertainty  is  not  known,  but  in 
practice  is  assumed  to  be  adequately  approximated  by  some  probability 
density  function  which  is  tractable  to  mathematical  analysis.  In  the  Monte 
Carlo  simulation  of  engagements  by  the  current  MGEM,  uncertainties  Nos. 
1-5  are  assumed  to  be  described  by  independent  Gaussian  distributions.  We 
include  uncertainty  No.  6  in  the  tabulation  even  though  the  firing  time  is 
assumed  to  be  precisely  executed  in  the  present  version  of  MGEM.  Such  a 
causal  (or  unitary)  PDF  is  described  in  the  mathematical  formulation  by  a 
Dirac  delta  function  (Reference  2).  Uncertainty  No.  7  is  usually  taken  into 
account  in  end-game  methodologies  by  using  prediction  algorithms  based 
on  noisy  estimates  of  the  target  state.  These  latter  result  from  uncertainty 
No.  8,  and  are  generated  in  the  code  by  imposing  appropriate  sensor  errors 
on  a  predetermined  precisely  known  flight  path  which  is  representative  of  a 
selected  aircraft  maneuver. 


4L.  Kruse,  and  P.  Brizzolara,  “An  Analytical  Model  for  Deriving  Conditional 
Probabilities  of  Kill  for  Target  Components,”  U.S.  Army  Ballistic  Research  Laboratories 
Report  No.  1563,  December,  1971,  (UNCLASSIFIED). 
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We  will  now  outline  a  brief  summary  of  a  typical  MGEM  calculation 
where  the  firing  doctrine  for  a  total  engagement  with  a  target  aircraft  is 
simulated.  This  outline  is  non-mathematical  and  is  used  to  introduce  the 
random  variables  in  the  context  of  which  they  will  be  used  later  to  con¬ 
struct  the  expected-value-of-kill  integrals  for  a  selected  set  of  conditions  in 
an  engagement. 

1.  Select  the  type  of  aircraft  maneuver  for  which  the  expected  value  of 
kill  is  desired,  and  retrieve  the  predetermined  flight  path  for  this  type  of 
maneuver. 

2.  Simulate  radar  fixes  on  the  target  aircraft  at  selected  times.  The 
location,  velocity,  and  acceleration  of  the  target  are  either  measured  or  cal¬ 
culated  for  these  times. 

3.  Using  the  data  derived  in  Step  2,  calculate  an  optimal  prediction  of 
future  target  locations.  Based  on  the  chosen  firing  doctrine,  select  a  value 
for  the  firing  time  T0  of  the  first  bullet  in  the  engagement. 

4.  The  location  r*,  of  a  bullet  at  any  time  t  is  assumed  to  be  adequately 
approximated  by  the  bullet  flyout  equations  given  in  the  next  section  of 
this  report.  Using  the  data  derived  in  Step  3  and  the  mean  values  of  Vm,  C4, 
Wr  Wv,  and  T0,  determine  an  estimated  interception  point  (r6)e  and  inter¬ 
ception  time  te  for  the  target  aircraft  by  the  bullet.  Values  for  the  polar 
aiming  angles  (©o^o)  are  also  derived  in  this  calculation. 

5.  Pick  a  set  of  values  for  the  initial  direction  (0o,*o)  °f  the  bullet  as  it 
leaves  the  gun  by  sampling  the  PDF  which  describes  the  uncertainties  asso¬ 
ciated  with  implementing  the  aiming  angles  calculated  in  Step  4.  In 
MGEM,  this  PDF  is  taken  to  be  Gaussian  where  the  values  of  the  aiming 
angles  (©o.’I'o)  are  taken  to  be  the  mean  values  of  (©o>*o)-  Additionally,  an 
error  bias  is  added  which  systematically  increases  the  aiming  angle  uncer¬ 
tainty  for  later  bullets  in  a  burst.  These  sample  values  are  identified  as 
[(^o).y.(^o).;]  where  i  is  the  running  index  over  bullets  in  an  engagement  and  j 
is  the  running  index  over  engagements.  The  number  of  bullets  fired  in  an 
engagement  is  taken  to  be  /  and  the  number  of  engagement  replications  in 
the  Monte  Carlo  calculation  is  taken  to  be  J. 

6.  Pick  a  value  for  each  of  the  random  variables  ( Vm, C4, by  sam¬ 
pling  the  appropriate  Gaussian  PDF.  These  sample  values  are  identified  as 
[( ed)ir( w*)iy( %)•;! •  These  values  are  used  below  to  calculate  the  necessary 
points  on  the  trajectory  of  the  sample  bullet. 


-12- 


7.  Simulate  the  firing  of  the  sample  bullet  and  calculate  its  true  loca¬ 
tion  (rfc)e  at  the  predicted  time  of  interception  te.  The  location  of  the  bullet 
which  has  been  in  flight  for  a  time  ( fy)e  is  calculated  by  applying  a 
differential  correction  for  the  effect  of  gravity  and  wind  to  the  line  of  sight 
(LOS)  range  vector  {(0o)»y>(^oM  (Figures  1  and  2). 

8.  Determine  the  true  location  (r,)4  of  a  reference  point  in  the  target 
aircraft  at  time  te  and  calculate  the  distance  between  (r,)e  and  (r6)e. 

9.  Calculate  the  distance  between  the  bullet  and  the  target  at  intervals 
of  At  about  te  until  a  minimum  distance  is  obtained. 

10.  Using  the  two  calculated  distances  adjacent  to  the  calculated 
minimum  distance,  interpolate  to  obtain  an  approximation  of  the  distance 
of  closest  approach  (miss  distance)  between  the  target  and  the  bullet.  The 
location  of  the  bullet  at  its  point  of  closest  approach  is  identified  as  (r6)m 
and  the  associated  time  is  identified  as  fTO. 

11.  The  orientation  of  the  target  aircraft  in  the  GCS  can  be  fixed  by 
either  locating  three  points  (the  aforementioned  reference  point  plus  two 
additional  points  in  the  aircraft),  or  by  two  vectors  emanating  from  the 
reference  point  (Figure  2).  Without  further  specification,  identify  these  two 
vectors  as  ua.  Determine  the  orientation  (u>,)m  of  the  target  at  time  tm.  Cal¬ 
culate  the  area  of  the  target  presented  to  the  incoming  bullet.  Represent¬ 
ing  this  presented  area  as  a  circle,  determine  if  the  bullet  impacts  the  tar¬ 
get.  Go  to  Step  14  if  the  bullet  misses  the  target. 

12.  Determine  if  the  aircraft  is  killed  by  comparing  a  random  number 
with  the  ratio  of  the  vulnerability  area  to  the  presented  area  of  the  aircraft. 

13.  Provision  is  also  made  to  treat  multiply-vulnerable  components, 
the  destruction  of  any  one  of  which  will  not  by  itself  result  in  target  kill.  If 
single-hit  kill  is  not  found  in  Step  12,  determine  by  random  number  draw 
whether  one  of  these  components  was  destroyed  and  whether  this  in  con¬ 
junction  with  previous  damage  results  in  target  kill. 

14.  Reiterate  Steps  1-13  for  the  remaining  bullets  in  the  burst  as  dic¬ 
tated  by  the  firing  doctrine.  If  target  kill  is  achieved,  terminate  the 
engagement.  Otherwise,  reiterate  burst  firing  until  the  firing  doctrine  dic¬ 
tates  engagement  termination. 
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Figure  2.  The  Vectors  (Gun-Coordinate  System)  Associated  with  an  Air-Defense 
End-Game  Problem 


15.  Reiterate  Steps  1-14  until  a  total  of  J  replications  of  the  engage¬ 
ment  have  been  obtained. 

16.  Calculate  the  expected  value  of  kill  E{K)  for  the  engagement  as  the 
average  kill  probability  of  the  J  replications,  that  is 

pttf\ _ Number  of  Kills  /.\ 

Number  of  Replications  ' 

Step  16  completes  the  outline  of  a  typical  Monte  Carlo  calculation  by 
MGEM.  We  have  discussed  it  in  some  detail  to  provide  a  concrete  example 
of  Monte  Carlo  methodology,  before  turning  to  more  general  mathematical 
formulations.  In  the  next  section  of  this  report,  we  will  first  use  the  uncer¬ 
tainties  tabulated  earlier  in  this  section  to  construct  the  expected-value-of- 
kill  integral  for  a  single  critical  component  for  a  single  set  of  firing  condi¬ 
tions  in  an  engagement.  We  will  then  restructure  the  integrand  in  this 
integral  by  changing  the  variables  of  integration  (0o»^o<ci)  to  a  new  set  which 
permits  the  picking  of  sample  bullets  which  will  with  a  higher  probability 
interact  with  the  target  aircraft.  Then,  in  Section  III,  we  will  outline  in 
turn  Monte  Carlo  procedures  which  can  be  used  to  calculate  a  forward  esti¬ 
mate  of  the  original  “natural”  expected-value-of-kill  integral  and  an 
adjoint,  or  backward,  estimate  of  its  restructured  form.  The  term  “for¬ 
ward  estimate”  is  used  here  to  describe  a  calculation  where  each  sample 
bullet  is  picked  at  the  gun  and  then  transported  toward  the  target  while 
the  term  “adjoint  or  backward  estimate”  is  used  to  describe  a  calculation 
where  each  sample  bullet  is  picked  at  a  point  on  its  trajectory  near  its  tar¬ 
get  and  then  transported  back  to  the  gun.  Finally,  in  Section  IV,  we  devise 
an  example,  two-dimensional  problem  and  outline  its  solution  by  using,  in 
turn,  forward  and  adjoint  Monte  Carlo  sampling.  The  calculational 
efficiency  of  forward  vs.  adjoint  estimation  is  compared  for  different  target 
sizes.  A  simplified  outline  of  the  Monte  Carlo  evaluation  of  expected-value 
integrals  is  given  in  the  appendix  at  the  end  of  this  report,  and  should  be 
read  before  proceeding  to  the  next  section. 
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II.  THE  CONSTRUCTION  OF  EXPECTED-VALUE-OF-KILL  INTEGRALS 


Our  main  interest  in  this  study  is  to  pick  sample  bullets  for  Monte 
Carlo  simulation  calculations  which  have  a  high  probability  of  interacting 
with  the  target  aircraft.  This  interaction  for  a  PD  bullet  will  require  a 
direct  hit  on  the  target,  while  a  possible  interaction  for  a  PX  bullet  will 
merely  require  that  the  bullet  pass  into  the  fuzing  region  which  surrounds 
the  target.  We  will  first  formulate  a  “natural”  form  of  the  expected-value- 
of-kill  integral  and  then  change  three  of  its  variables  of  integration  so  that 
the  transformed  (adjoint)  integral  is  in  a  form  which  is  amenable  to  picking 
sample  bullets  which  will  interact  with  a  high  probability  with  the  target. 

Since  our  principal  interest  here  is  the  development  of  alternate  formu¬ 
lations  of  the  expected-value-of-kill  integrals  and  a  comparison  of  the 
expected  efficiencies  of  the  Monte  Carlo  evaluation  of  these  different  formu¬ 
lations,  we  will  not  describe  the  detail  associated  with  locating  a  target  and 
firing  a  bullet  to  obtain  its  kill,  but  will  merely  identify  the  various  associ¬ 
ated  uncertainties  and  represent  them  by  probability  density  functions 
(PDF’s).  The  general  form  of  the  expected  value  of  kill  E{K }  of  a  critical 
component  in  a  target  aircraft  by  a  bullet  fired  from  a  gun  at  time  /0  is 
given  by  the  integral  of  the  product  of  the  various  PDF’s  used  to  predict 
the  uncertainties  associated  with  locating  the  target,  aiming  the  gun,  firing 
the  bullet,  transporting  the  bullet  toward  its  aircraft  target,  and  a  kill 
function  which  approximates  the  effect  of  a  bullet  on  the  operation  of  a 
critical  component  in  the  aircraft.  In  order  to  provide  a  similar  form  for 
both  the  natural  and  adjoint  formulation  of  both  PD  and  PX  munitions, 
we  define  a  region  /?f  which  for  a  PD  bullet  barely  encloses  the  volume  of 
the  target  and  for  a  PX  bullet  barely  encloses  the  volume  of  the  target  plus 
the  associated  region  about  the  target  in  which  a  fuzing  action  by  the  bul¬ 
let  can  occur.  We  will  choose  a  sphere  which  seems  to  be  the  most  con¬ 
venient  form  for  calculations  (Figure  3).  In  each  case,  an  associated  surface 
Sf  is  identified  which  is  that  part  of  the  boundary  of  /?f  which  can  be  inter¬ 
sected  by  bullet  trajectories. 

We  assume  that  the  lumping  of  the  various  uncertainties  associated 
with  an  air-defense  system  into  the  form  described  earlier  in  Section  I  of 
this  report  will  permit  the  expected-value-of-kill  integral  to  be  formulated 
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E{K{  r0,e0,v0,e0,v0,  vm,cit  vvF)}  = 

/ *RV /  '4(<oWo.^o)ct^o.V»ol^o.V»o)  Sin0o  E{vm)  G(ei)  U{w*’wv)  X 

QfGI(I‘«)r(fa)p(w«)r(*,*)p('*)r(£J*)i]  dw x  dwv  dvm  dcd  d0O  #0  d0 0  #0  ^  •  (2) 

The  quantity  R%  is  the  region  in  (to,0o,ii>Q,OQ,il>o,vm,c4,ivtywi)  space  which 
includes  all  possible  firing  events  where  the  bullet  trajectories  will  intersect 
Sv  In  practice,  since  the  commonly  used  Gaussian  PDF’s  are  truncated  so 
as  to  include  only  a  reasonable  region  about  their  mean  value,  we  will  take 
/?$  to  include  only  those  points  where  all  of  the  various  PDF’s  are  non-zero. 

This  integral  can  be  evaluated  when  the  physics  describing  the  bullet 
flyout  is  sufficiently  well  understood  so  that  a  description  of  all  bullet- 
trajectory  intersections  of  Sf  is  calculable  to  an  adequate  accuracy.  For 
example,  in  an  end  game  methodology  developed  by  Ford5,  the  bullet  flyout 
to  target  impact  is  calculated  by  using  iterative  procedures  to  solve  the 
differential  equation  constructed  from  the  forces  acting  on  a  bullet  in  flight. 
In  the  present  study,  after  defining  the  quantities  used  above,  we  will  apply 
the  analytic  functions  used  by  Reference  1  to  predict  the  flyout  of  bullets. 

The  order  of  integration  of  the  integrals  used  in  this  study  is  assumed 
to  start  at  the  innermost  integral  and  proceed  toward  the  outer  integrals. 
In  general,  we  express  the  arguments  of  functions  in  terms  of  their  most 
convenient  variable  association  without  regard  to  whether  the  variables 
used  are  variables  of  integration.  Needless  to  say,  before  conducting  an 
integration,  all  functions  must  be  given  in  terms  of  the  variables  of  integra¬ 
tion. 

In  Equation  2,  we  have  allocated  the  representation  of  the  various 
uncertainties  to  six  probability  density  functions  where  each  function 
describes  a  particular  aspect  of  the  total  problem.  The  uncertainties  associ¬ 
ated  with  the  firing  time  are  described  by  the  PDF  A(t0).  If,  as  presently 
assumed  in  MGEM,  a  precisely  known  firing  time  is  assumed,  its  PDF  is  a 
Dirac  delta  function,  that  is 

(3A) 


&“DIVAD  Contract  Data  Requirement  No.  A00101,  System  Simulation  Analyst  Manual, 
Revision  A,”  Ford  Aerospace  and  Communication  Corporation,  Aeronutronic  Division, 
Newport  Beach,  California  92663,  November  1980,  (Competition  Sensitive). 
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where  fc  is  the  precisely  known  firing  time.  The  remaining  temporal  vari¬ 
ables,  which  will  be  explicitly  used  in  the  adjoint  formulation  later  in  this 
section,  are  related  by 

<  =  <o  +  </•  (3B) 

where  tf  is  the  time  of  flight  of  a  bullet  and  t  is  the  target  time.  For  a 
given  value  of  T0,  the  associated  infinitesimal  times  are  related  by 

dtj  =  dt.  (3C) 

The  uncertainties  associated  with  locating  the  target,  predicting  its 
future  locations,  and  calculating  the  aiming  angles  for  the  gun  are  described 
by  the  joint  PDF  B(0o,rpo).  In  practice,  the  calculated  values  (tfo.V’o)  of  the 
aiming  angles  are  derived  from  an  optimal  estimate  of  future  target  loca¬ 
tions  based  on  a  sequence  of  radar  fixes  obtained  prior  to  the  firing  of  the 
bullet.  An  uncertainty  is  associated  with  each  radar  fix  so  that  an  estimate 
of  future  locations  of  the  same  target  derived  individually  from  two  alter¬ 
nate  sequences  of  radar  fixes  will  not  in  general  be  the  same.  Consequently, 
(0O,V>O)  is  generally  some  function  of  the  radar  uncertainties  and  the  target 
maneuver  whose  analytic  form  would  probably  be  very  difficult  to  derive 
and  very  cumbersome  to  write. 

In  addition  to  the  aiming  uncertainties  associated  with  (0O,^>O),  the  accu¬ 
racy  of  the  true  firing  direction  (0O,^O)  °f  bullets  is  further  degraded  by 
uncertainties  associated  with  positioning  the  gun  for  firing.  We  assume 
here  that  the  conditional  probability  density  C[e0,xp0 |?0,^0)  of  real  firing 
angles  is  the  probability  density  per  unit  solid  angle  Cl.  If  the  air-defense 
system  has  a  non-zero  probability  for  firing  into  the  upper  hemisphere 
(atmosphere)  and  is  zero  elsewhere,  then  its  normalization  is  given  by 

W 

2*2  ^ 

/  C^o^oK^o)^  =  J  /  Woi'PoWoi'Po)  Sinffo  Mo  =  1.  (3D) 

2  if  0  0 

for  any  given  set  of  values  for  (?0,^0).  For  example,  a  constant  value  of  l/2n 
for  Cf0o>V#oA’o)  would  correspond  to  firing  isotropically  into  the  upper  hem¬ 
isphere. 

The  variation  in  muzzle  velocities  among  the  bullets  is  described  by  the 
PDF  f\vm).  The  value  for  the  drag  coefficient  C4  is  taken  to  be  constant 
during  the  flyout  of  a  bullet,  and  its  values  for  an  aggregate  of  bullets  is 
predicted  by  the  PDF  G(ci).  The  values  for  the  components  of  force  due  to 
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wind  are  predicted  by  the  joint  PDF  U[w„ wy)  where  V  in  MGEM1  is 
assumed  to  be  separable  into  the  two  functions 

U(w,,wv)  =  /,( wt)  f2(wt).  (3E) 

A  set  of  values  when  picked  for  Wg  and  W9  is  also  taken  to  remain  constant 
for  each  bullet  during  its  flyout. 

As  noted  earlier4,  the  loss  of  capability  by  a  target  component  due  to 
the  interaction  of  a  bullet  is  generally  approximated  as  a  kill-probability 
function  of  certain  random  variables  which  can  be  used  to  characterize  an 
aggregate  of  bullets  or  metal  fragments  as  they  impact  the  outer  surface  of 
a  critical  component.  The  formulation  used  in  Equation  2  differs  from  the 
definition  given  by  Reference  4  in  that  the  arguments  of 
^?[(r«)p(^«)r(w,)P(r*)r(r*)r(wt)<l  are  the  values  in  the  GCS  of  the  random  vari¬ 
ables  at  the  time  when  a  bullet  or  fragment  impacts  a  hypothetical  surface 
S(  instead  of  the  real  exterior  surface  of  a  target  component.  As  noted  ear¬ 
lier  in  Section  1,  the  position  vectors  r4  and  rb  are  used  to  locate  a  reference 
point  in  the  target  aircraft  and  bullet,  respectively,  and  the  vectors  r,  and 
Tb  are  the  velocities  of  these  reference  points  in  the  target  and  bullet, 
respectively.  The  vector  pairs,  u>,  and  ub,  define  the  orientation  of  the  bul¬ 
let  and  target,  respectively,  where  all  vectors  are  defined  in  the  GCS.  The 
subscript  f  is  used  to  identify  the  associated  vector  at  the  intersection  of 
the  bullet  with  Sf. 

The  quantity  Qf  is  itself  the  expected  value  of  kill  of  the  bullet  at  its 
intersection  of  Sf,  and,  in  principle,  assuming  that  the  location  and  velocity 
vectors  provide  an  adequate  characterization  in  a  stochastic  formulation  of 
bullet  flyout,  could  have  been  described  in  Equation  2  by  increasing  the 
dimensionality  of  the  multiple  integrals  to  include  the  appropriate  functions 
which  describe  the  further  transport  of  the  bullet  into  f  and  its  interaction 
with  the  target  aircraft.  In  practice  during  current  Monte  Carlo  calcula¬ 
tions  for  PD  bullets,  the  kill  function  is  approximated  at  the  surface  of  the 
aircraft  for  damage  due  to  both  the  penetration  of  detonation  fragments 
into  the  target  and  the  transmission  of  shock  to  components  distant  from 
the  point  of  impact.  In  practice  for  a  PX  bullet  in  forward  Monte  Carlo 
calculations,  an  estimate  of  Qf  would  be  obtained  for  each  firing  event  by 
transporting  the  bullet  into  Rf,  simulating  its  detonation  or  non-detonation, 
and  then  deriving  the  kill  probability  from  a  kill  function  such  as  that 
described  in  Reference  4. 
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We  will  assume  here  that  the  mathematical  relations  used  by  Reference 
1  can  be  used  to  adequately  predict  the  flight  of  bullets  toward  a  target  air¬ 
craft.  According  to  Reference  1,  McShane,  Kelley,  and  Reno®  have  derived 
the  LOS  range  /  to  be  given  by  the  Siacci  formula, 

l  =  Vj{ 

1  +  ca/^  l! 

where  is  the  time  that  a  bullet  has  been  in  flight.  The  effect  of  wind  and 
gravity  are  assumed  to  be  differential  corrections  to  the  predicted  LOS 
range  vector  (O0,ip0,f)  which  are  given  in  Cartesian  coordinates  by: 

Ax  =  wt  tj  - 

Ay  =  u-„  tj- 

A  _  »2  Z+'i\fc  */ 

The  quantities,  Ax,  Ay,  and  As,  are  the  differential  corrections,  and  g  is  the 
acceleration  due  to  gravity.  ~~ 

As  outlined  in  the  next  section  of  this  report,  a  forward,  or  “natural”, 
Monte  Carlo  evaluation  of  the  multiple  integrals  in  Equation  2  can  be  con¬ 
ducted  by  picking  values  for  each  of  the  random  variables  of  integration  by 
sampling  the  appropriate  PDF’s  to  obtain  sample  bullet  trajectories.  Each 
sample  bullet  is  then  transported  along  its  trajectory  using  the  preceding 
relations  to  determine  if  it  interacts  with  the  target  aircraft.  The  effect  of 
each  interaction  is  determined  and  the  expected  value  of  kill  is  calculated 
as  the  mean  of  a  large  number  of  similar  events.  However,  using  calcula- 
tional  procedures  such  as  these  for  cases  where  the  target  is  at  a  large  dis¬ 
tance  from  the  gun,  a  large  fraction  of  the  bullets  will  miss  the  target  and 
contribute  only  zeroes  to  the  expected  value  of  kill  probability.  Conse¬ 
quently,  the  efficiency  of  the  forward  mode  decreases  as  the  hit  probability 
decreases  until,  for  very  low  hit  probabilities,  a  prohibitive  number  of  sam¬ 
ple  histories  would  be  required  to  obtain  an  acceptable  fractional 

®E.J.  McShane,  J.L.  Kelley,  and  F.V.  Reno,  Exterior  Ballistics.  University  of  Denver 
Press,  1953. 
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uncertainty  — =—  (where  SPH  is  the  standard  deviation  of  a  hit-probability 

f*H 

estimate  PH). 


We  will  now  attempt  to  transform  the  “natural”  expected-value-of-kill 
integral  (Equation  2)  to  a  form  such  that  only  bullets  which  have  a  high 
probability  of  interacting  with  the  target  can  be  picked.  As  a  first  step,  we 
note  that  the  region  R t  is  the  smallest  sphere  which  just  encloses  the  target 
aircraft  when  the  bullets  are  PD  and  which  just  encloses  the  target  aircraft 
plus  its  associated  fuzing  region  when  the  bullets  a_re  PX.  The  radius  of  the 
sphere  in  each  case  is  taken  to  be  the  constant  a.  We  construct  a  Target 
Coordinate  System  (TCS)  (Figure  4)  whose  origin  is  at  the  center  of  R(  and 
whose  orientation  is  such  that  the  Cartesian  coordinates  of  each  point  in 
the  TCS  is  related  to  the  same  point  in  the  GCS  by  a  simple,  time- 
dependent  translation  of  the  center  (*a,ya,2a)  of  R(  in  the  GCS;  that  is,  the 
( x ’  ,y '  ,z’  )  axes  in  the  TCS  are  parallel  to  and  have  the  same  orientation 
as  the  (x,i i,z)  axes  in  the  GCS  (Figure  5).  Similarly,  the  spherical  coordi¬ 
nates  in  the  TCS  are  similar  to  and  have  the  same  orientation  as  the 

spherical  coordinates  (0,^,r)  which  are  used  in  the  GCS.  The  surface 
which  bounds  /?f  is  in  the  TCS  just  that  part  of  a  —  a  which  can  be 
impacted  by  bullets  fired  at  the  origin  of  the  GCS. 

We  next  note  that  a  large  fraction  of  the  bullet  trajectories  which  inter¬ 
sect  Sf  will  interact  with  the  target  aircraft;  in  fact,  this  fraction  is  rela¬ 
tively  insensitive  to  factors  affecting  the  hit  probability  such  as  target  size 
and  range.  If  we  could  change  some  of  the  variables  of  integration  used  in 
Equation  2  to  a  new  set  of  variables  such  that  all  sample  trajectories  will 
intersect  Sf,  then  we  should  expect  that  more  efficient  calculations  of  esti¬ 
mates  than  are  currently  being  made  could  be  developed  for  air-defense 
problems  where  the  expected  value  of  hit  is  small. 

To  accomplish  this  change  of  variables,  we  note  that  the  location 
[(£»)(>( ViM**);]  in  the  GCS  of  the  bullet  intersection  of  Sf  is  related  to  its  loca¬ 
tion  [( V  )r(y*'  )r(V  y  in  the  TCS  by: 


(**) f  =  (*b  )<  +  *.(0. 

(5A) 

(y*)f  =  (Vb  )f  +  va(<), 

(5B) 

(*»). f  =  (V  )f  +  *«(0, 

(5C) 

where  (*.(<), ya( /).*•(  01  is  the  time-dependent  location  of  the  center  of  /?f  in  the 
GCS.  These  relations  become  in  spherical  coordinates: 
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gure  4.  The  Target  Spherical  Coordinate  System  Associated  with  the 

Change  of  Integration  Variables 
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SkkLL  v. 


l(  sin  O0  cos  ip o  +  WJ 
/(  sin  0o  sin  V’o  +  tuyl 


«/)<  - 


t'A-v- 


/f  cos  0O  - 


*'A2 

•  * 

3  +  Cd(tf)( 

6 

1  + 

=  a  sin^i  cos{fc  +  za(t), 
=  a  sin<^  sinf6  +  yjt), 

=  a  cos <f>b  +  *„(<). 


(6A) 

(6B) 

(6C) 


The  quantities  /f  and  (fy)f  are  the  LOS  range  and  the  bullet  time-of-flight, 
respectively,  for  trajectory  intersections  of  Sr  Further  noting  that  /f  is  given 
by  the  Siacci  equation, 

l’m  (</)f 


1  +  cd  (*/)f 


(4A) 


and  the  temporal  variables  are  related  by, 


t  —  t0  +  (ty)f ,  (3B) 

Equations  6A-6C  provide  three  independent  equations  relating  three  of  the 
original  variables  of  integration,  namely  the  spatial  variables  (fy^o)  ar>d 
either  vm  or  cd,  and  the  TCS  spatial  variables  (<*,£)  and  the  temporal  vari¬ 
able  t. 


We  will  arbitrarily  select  cd  to  participate  in  the  transformation  so  that 
the  sought  change  of  variables  is  from  (0o,tpo,cd)  to  (<>,£,0  Assuming  that  all 
points  in  R%  will  transform  one-to-one  to  a  region  in 
space,  we  can  write  the  transformed  expected-value-of-kill  integral  as 

E{ K(  ro,0o,¥ 0-e0,*o.  Em.Q,  W„  W,)}  = 

E{/q  T0,©0,^0,4>,H,  r,  Vm,  W„  W,]}  = 


/•  *T*J  ^(<0)  B(®o.  V’o)  C^o.^ofo.^o)  Sin  °o  E[vm)  G(cd)  U[wt,wf)  x 

Q(rKr»)r('*)r(u,«)fl  dw •  dwt  dv  dv*  #o  <*<o  •  (7A) 

In  this  formulation,  we  have  assumed  that  the  component  kill  function  Q* 
is  given  in  the  TCS.  The  quantity  T  is  taken  to  be  the  random  variable 
whose  values  are  t.  The  Jacobian  J£<t>,Z,t)  is  the  absolute  value  of  the  deter¬ 
minant  of  the  indicated  partial  derivatives  at  the  time  when  the  bullet  tra¬ 
jectory  intersects  Sf,  that  is 
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d0o  d0o  dOQ 

~d$  Hz  HT 

dV'o  dip  Q  dtpQ 
dtp  HT  dt 

dcj  dcj  dcj 

HT  Hz  Hi 

As  is  always  the  case,  all  functions  in  the  integrand  must  be  given  in  terms 
of  the  variables  of  integration  before  conducting  an  integration. 

Equation  7A  is  in  a  form  where  our  sampling  objectives  can  be 
attained.  A  point  on  S(  at  some  particular  target  time  in  conjunction  with 
the  known  location  of  the  gun  and  a  set  of  sample  values  for  the  remaining 
variables  of  integration,  unchanged  from  the  forward  formulation  (Equation 
2  ),  will  completely  determine  a  bullet  trajectory.  A  Monte  Carlo  evalua¬ 
tion  will  be  developed  in  the  next  section  of  this  report  where  sample  tra¬ 
jectories  are  picked  using  these  changed  variables  of  integration. 

However,  we  will  note  here  that  the  calculational  efficiency  may  be 
improved  by  picking  the  sample  points  (r4)f  directly  from  S(  instead  of  first 
picking  values  for  the  angles  (4>,S)  and  then  locating  the  associated  point  on 
Sf.  The  relevant  quantities  are  related  by 

dS<  =  a  sin<£  d<p  dZ .  (8) 

Procedures  for  conducting  such  sampling  as  well  as  an  analysis  of  their  cal¬ 
culational  efficiency  should  be  investigated.  However,  we  will  not  expand 
further  on  this  line  of  development  in  the  present  study. 


WM  = 


d(0o,tl>o,cj) 

d(4,Z,t) 
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HI.  A  GENERAL  MONTE  CARLO  EVALUATION  OF  THE  NATURAL 
AND  THE  TRANSFORMED  EXPECTED-VALUE-OF-KILL  INTEGRAL 


A  Monte  Carlo  procedure  for  evaluating  the  multiple  integrals  in  equa¬ 
tion  2  is  outlined  below  and  illustrated  in  Figure  6.  In  order  to  keep  the  fol¬ 
lowing  outline  as  general  as  possible,  we  will  neither  describe  the  sampling 
procedures  used  to  pick  variable  values  from  the  various  PDF’s  nor 
attempt  to  develop  optimum  procedures  for  picking  values  for  the  new 
variables  R,<^,£).  This  procedure  is  composed  of  the  following  steps: 

1.  Select  the  type  of  aircraft  maneuver  for  which  the  expected  value  of 
kill  is  desired,  and  retrieve  the  predetermined  flight  path  for  this  type  of 
maneuver. 

2.  Simulate  radar  fixes  on  the  target  aircraft  at  selected  times.  The 
location,  velocity,  and  acceleration  of  the  target  are  either  measured  or  cal¬ 
culated  for  these  times. 

3.  Use  the  data  just  derived  in  Step  2  to  calculate  an  optimal  predic¬ 
tion  of  future  target  locations.  Applying  some  preselected  firing  doctrine, 
pick  a  value  for  the  firing  time  T0  by  sampling  This  sample  firing  time 
is  identified  as  (/0),. 

4.  The  location  rb  of  a  bullet  at  time  t  is  assumed  to  be  adequately 
approximated  by  the  flyout  relations  (Equations  4A-4D)  given  in  Section  II 
of  this  report.  Use  the  data  derived  in  Step  3  and  the  mean  values  of 
(VmyCd,WrWy)  to  determine  an  estimated  interception  point  (rt)e,  intercep¬ 
tion  time  t(,  and  aiming  angles  [(?0),»(^o)J-  The  reader  should  note  that  this 
picking  of  a  set  of  sample  values  for  the  aiming  angles  is  accomplished 
without  ever  knowing  the  explicit  form  of  the  PDF  fl(0o,V>o). 

5.  Pick  a  set  of  values  for  the  initial  direction  (©o-^o)  bullet  as  it 

leaves  the  gun  by  sampling  C}0o,V>o](^o)..(V;o)i}  sin  •  These  sample  values  are 
identified  as  ((0o).-,(^o)J-  Procedures  for  picking  sample  values  for  a  random 
variable  from  its  PDF  are  discussed  by  Shreider.7  These  sample  values  plus 
the  fact  that  the  birth  site  of  each  bullet  is  located  at  the  origin  of  the  GCS 


7N.P.  Buslenko,  D.I.  Golenko,  Yu.A.  Shreider,  I.M.  Sobol’,  and  V.G.  Sragovicb,  The 
Monte  Carlo  Method.  The  Method  of  Statistical_Trials.  Yu.A.  Shreider,  Editor, 
Translated  from  the  Russian  by  G.J.  Tee,  Translation  edited  by  DM.  Parkyn, 
Pergamon  Press,  Aylesbury,  Bucks,  Great  Britain,  1967. 
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Figure  6.  The  Forward,  or  “Natural”,  Monte  Carlo  Estimation  of  the 
Expected- Value-of-Kill  of  a  Critical  Component  in  a  Target 
Aircraft 
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provides  the  initial  conditions  for  the  bullet  flyout. 

6.  Pick  a  value  for  each  of  the  random  variables  (Vm,Cd,W„Wy)  by  sam¬ 
pling  the  appropriate  PDF.  These  sample  values  are  identified  as 
l(Oii(cA«W»»(®»)d-  These  values  are  used  below  to  calculate  the  necessary 
points  on  the  trajectory  of  the  sample  bullet. 

7.  Simulate  the  firing  of  the  sample  bullet  and  calculate  its  true  loca¬ 
tion  (rt)t  at  the  predicted  time  of  interception  te.  The  location  of  a  bullet 
which  has  been  in  flight  for  a  time  (fy)e  is  calculated  by  applying  the 
differential  corrections  (Equations  4B-4D)  to  the  line  of  sight  (LOS)  range 
vector  l(0oMV>oWe]  (Equation  4A). 

8.  Determine  the  true  location  (r#)e  and  orientation  (u?a)e  (Figure  2)  of 
the  target  aircraft  at  time  te. 

9.  Reiterate  Step  8  at  selected  intervals  of  At  about  te  until  an  intersec¬ 
tion  with  S(  is  bracketed  (if  such  occurs).  Set  the  kill-probability  score  to  0 
and  go  to  Step  14  if  the  bullet  does  not  impact  S{. 

10.  Using  the  bracketing  data  acquired  in  Step  9,  interpolate  to  find  an 
adequate  approximation  of  the  bullet  intersection  of  Sr  The  location  of  the 
bullet  at  this  point  in  its  trajectory  is  identified  by  the  position  vector  (rt)f 
and  the  associated  time  is  t{.  The  bracketing  data  acquired  in  Step  9  can 
also  be  used  to  calculate  an  adequate  approximation  of  the  bullet  velocity 
at  time  t(  if  this  quantity  is  needed  for  the  calculation  of  kill  probabilities. 
Continue  to  Step  11  for  a  PD  bullet,  otherwise  go  to  Step  12  for  a  PX  bul¬ 
let. 


11.  Calculate  the  kill  probability  for  each  critical  component  and  then 
go  to  Step  13.  Each  kill-probability  score  is  identified  as  X,y  where  j  is  the 
running  index  over  critical  components  and  J  is  the  number  of  critical  com¬ 
ponents  in  the  target  aircraft. 

12.  Determine  the  time  when  a  PX  bullet  impacts  the  fuzing  region, 
and  start  transporting  the  bullet  into  this  region.  Track  the  bullet  until 
either  a  fuzing  action  occurs,  the  bullet  impacts  the  target  aircraft,  or  the 
bullet  leaves  the  fuzing  region.  The  burst  of  metal  fragments  produced  by  a 
detonation,  if  such  an  event  occurs,  must  then  be  transported  into  the  tar¬ 
get  and  the  kill  probabilities  of  impacted  critical  components  must  be 
determined.  A  detailed  description  of  such  calculations  is  beyond  the  intent 
of  this  study  and  will  not  be  addressed  here.  A  possible  approach  is  to 
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extend  the  stochastic  description  of  the  vulnerability  process  suggested  by 
Beverly8®  for  critical  components  exposed  to  spall  fragments  in  heavily 
armored  land  vehicles  (tanks)  to  include  the  higher  target  velocities  associ¬ 
ated  with  fixed-wing  aircraft.  We  have  not  attempted  to  specify  the  detail 
with  which  the  target  aircraft  is  described  but  have  merely  assumed  here 
and  in  Step  11  that  it  is  adequately  described. 

13.  Accumulate  each  score  in  a  bin  reserved  for  the  associated  com¬ 
ponent.  Square  each  score  and  accumulate  each  squared  score  in  another 
bin  reserved  for  the  associated  component. 

14.  Reiterate  Steps  1-13  until  the  histories  for  a  total  of  I  sample  bul¬ 
lets  have  been  simulated. 

15.  Applying  the  Law  of  Large  Numbers2,  calculate  the  expected  value 
of  kill  for  each  critical  component  as 

E[K)j  =4  E  =  V  (9A) 

1  i  -  1 


16.  Calculate  an  estimate  SE{K} i  of  the  standard  deviation  of  the 
estimated  expected  value  of  kill  for  each  component  as10 


6E{K}i  = 


E  x5/-;xy 


(9B) 


Applying  the  Central  Limit  Theorem2,  this  standard  deviation  can  be  inter¬ 
preted  as  implying  that  approximately  68  percent  of  a  large  number  of 
similar  expected-value-of-kill  estimates  will  fall  within  one  standard  devia¬ 
tion  of  the  associated  universe  expected  value. 


8W.B.  Beverly,  “A  Stochastic  Representation  of  the  Vulnerability  of  a  Critical 
Component  in  a  Military  Vehicle  to  Metal  Fragments,”  Journal  of  Ballistics,  Vol.  5,  No. 
4,  October,  1981. 

°W.B.  Beverly,  “The  Application  of  the  Monte  Carlo  Method  to  the  Solution  of  the 
Interna]  Point  Burst  Vehicle  Ballistic  Model,”  Ballistic  Research  Laboratory  T  .nical 
Report  No.  02353,  August  1981,  (UNCLASSIFIED). 


10Y.Beers,  Introduction  to  the  Theory  of  Error.  Addison-Wesley  Publishing  Company, 
Inc.,  Reading,  Mass.,  1962. 
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Figure  7.  The  Adjoint  Monte  Carlo  Estimation  of  the  Expected-Value-of- 
Kill  of  a  Critical  Component  in  a  Target  Aircraft 
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Step  16  completes  the  outline  of  a  forward  Monte  Carlo  solution  of  the 
expected-value-of-kill  integral  used  in  Equation  2.  The  alternate  adjoint  for¬ 
mulation  given  by  Equation  7A  is  in  a  form  such  that  sample  bullets  which 
usually  impact  Sf  can  be  easily  picked.  A  general  Monte  Carlo  procedure  for 
evaluating  the  adjoint  integral  is  composed  of  the  following  steps  (Figure 

7)- 


1.  Select  the  type  of  aircraft  maneuver  for  which  the  expected  value  of 
kill  is  desired,  and  retrieve  the  predetermined  flight  path  for  this  type  of 
maneuver. 

2.  Simulate  radar  fixes  on  the  target  aircraft  at  selected  times.  The 
location,  velocity,  and  acceleration  of  the  target  are  either  measured  or  cal¬ 
culated  at  each  of  these  times. 

3.  Use  the  data  derived  in  Step  2  to  calculate  optimal  predictions  of 

future  target  locations.  Applying  some  preselected  firing  doctrine,  pick  a 
value  for  the  firing  time  T0  by  sampling  This  firing  time  is  identified 

as  (f0),  where  «  is  the  sample  index,  and  /  is  the  number  of  sample  histories 
to  be  used  in  the  Monte  Carlo  calculation. 

4.  Using  the  data  derived  in  Step  3  and-  the  _mean  values  of 
(Vm,Cd.  Wr  H’,),  calculate  values  for  the  aiming  angles  (©0-^o)  These  calcu¬ 
lated  values  are  identified  as  [(?<,), ,(V’o)il 

5.  Pick  a  set  of  sample  values  [( vm)iy (te,)„( u*y)J  for  ( Um,  W„  H^)  by  sampling 
the  associated  PDF’s. 

6.  We  now  introduce  the  first  step  in  the  adjoint  evaluation  which 
differs  from  the  forward  evaluation  just  outlined.  The  range  of  t ,  and 
correspondingly  (ty)f,  over  which  the  major  fraction  of  the  bullets  will 
impact  Sf  must  be  estimated  by  using  the  data  derived  in  the  first  five 
steps.  If  Gaussian  distributions  are  used  to  predict  the  uncertainties,  non¬ 
zero  probability  tails  will  extend  out  to  infinite  values  so  that  small,  but 
non-zero,  probabilities  will  exist  for  implausible  firing  or  bullet-flyout  condi¬ 
tions.  In  practice,  these  tails  should  be  truncated  and  a  reasonable  target- 
time  range  (/?,).  be  determined  which  will  include  only  plausible  events.  The 
width  of  this  sampling  time  range  is  identified  at  (AT),,  and  will  be  used  as 
a  factor  in  the  calculation  of  the  event  score. 
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7.  Pick  a  sample  value  f,  for  the  target  time.  We  will  assume  here  that 
the  sample  value  is  picked  with  equal  probability  at  any  point  on  [Rt){. 
Using  the  notational  convention  adopted  earlier,  the  associated  time  of 
flight  to  a  bullet  impact  with  Sf  would  be  identified  as  However,  for 
the  sake  of  simplicity,  we  will  take  the  quantity  (fy),  in  the  remainder  of  the 
adjoint  procedure  to  be  the  time  of  flight  of  the  «**  sample  bullet  to  its 
impact  with  Sf. 

8.  Locate  the  target  aircraft  at  f,  and  determine  in  the  gun  coordinate 
system  its  velocity  r ,(<,)  and  orientation  w4(/,).  Usejthis  data  to  locate  the 
spherical  region  (/ty,  and  its  associated  center  [(*,)r(ya)r(*,)f],  in  the  GCS.  If 
practical,  determine  the  boundary  of  the  associated  as  that  part  of  the 
outer  boundary  of  which  can  be  impacted  by  any  plausible  bullet  tra¬ 
jectories  having  the  variable  sample  values  just  picked  in  Steps  1-8. 


9.  We  are  now  ready  to  pick  a  sample  set  of  values  for  the  target- 
coordinate-system  angles  ($,=).  If,  as  anticipated,  the  exact  boundary  of  (5),- 
is  difficult  to  derive  or  its  shape  introduces  sampling  difficulties,  then  the 
following  procedure  may  be  used: 

(a)  Determine  the  range  (fy),  of  ^-values  and  the  range  (R^){  of  ^-values 
where  these  ranges  include  all  possible  target  hits. 

(b)  Using  these  ranges,  calculate  the  associated  widths  (A^),  and  (Ac),. 
These  widths  will  be  used  as  factors  in  the  calculation  of  the  event  score. 

(c)  Pick  a  sample  value  for  $  with  equal  probability  at  any  point  on  (/?^),. 
Similarly,  pick  a  sample  value  for  5  with  equal  probability  at  any  point  on 

This  set  of  values  is  identified  as 

(d)  Convert  the  sample  angles  just  picked  to  the  associated  point 
[(V  )PW  MV  )<],•  (TCS)  on  the  surface  of  (Rf),-. 

(e)  Determine  if  this  point  is  located  at  the  entrance  or  exit  of  the  sample 
bullet  trajectory  into 

(f)  Set  the  score  to  0  and  go  to  Step  17  if  the  point  is  a  trajectory  exit. 
Proceed  to  Step  10  if  the  point  is  a  trajectory  entrance  into  By 
definition,  (S^),  includes  only  trajectory  entrance  points  into  (/?  ),-. 


10.  Calculate  the  associated  set  of  values  for  (©o^o^)  by  inserting  the 
Siacci  Equation 


(O. •(*;).• 

1  +  ei  \/(vm).  (*/). 


(4A) 


and  the  relation, 
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=  (to),  +  (tf), 

between  temporal  variables  into  Equations  6B-6D: 


/  sin0n 


COS +  («>,),  (/;),  -  ~y-J  =  [(V  )f],  +  *.(*,)  , 

/sin0o  sin 0o  +  («',),  (</),■  -  ~yj  =  [(y*'  )f),  +  fjft)  , 


/  COS<?0  - 


Q  (*/)?  3  +  cd  (//), 


6  1  +  «<(</).  V/KX 


=  UV  )(),-  +  *.(<•■) 


(3B) 


(6A) 


(6B) 


(6C) 


The  resulting  set  of  simultaneous  equations  can  be  solved  by  using  an  itera¬ 
tive  procedure  such  as  Newton’s  method.11  This  set  of  values  is  identified 
as  [(fyMV’oM  <■*),]• 

11.  Calculate  the  killing  probability  of  the  sample  event  as  described 
eaHier  in  this  Section  in  the  outline  of  the  forward  Monte  Carlo  solution.  If 
the  bullet  velocity  is  needed  to  perform  this  calculation,  then  an  adequate 
approximation  can  be  obtained  from  the  change  in  the  bullet  coordinates 
for  a  small  increase  (or  decrease)  in  the  time  of  flight  ( t/){.  This  quantity  for 
either  case  is  identified  as  Q{,  and  will  be  used  as  a  factor  in  the  calculation 
of  the  event  score. 

12.  Calculate  the  quantity  Ct(^o)(,(0o)If(?o)1,(0o)1]  S«n(0o).-.  This  quantity  is 
identified  as  C,sin(0o)„  and  will  be  used  as  a  factor  in  the  calculation  of  the 
event  score. 

13.  Calculate  the  value  of  the  Jacobian  4(0,-,  {„(*/),]  (Equation  7B)  for 
the  event  where  the  necessary  derivatives  are  derived  from  Equations  6A- 
6C.  The  time  dependence  of  r,  is  derived  for  the  particular  target 
maneuver  used.  This  quantity  is  identified  as  J{,  and  will  be  used  as  a  factor 
in  the  calculation  of  the  event  score. 

14.  Calculate  the  value  of  the  PDF  Cf(c,),]  for  the  event.  This  quantity 
is  identified  as  G{,  and  will  be  used  as  a  factor  in  the  calculation  of  the 
event  score. 


^Peter  Henrici, 
N.Y.,  1664. 


iiS,  John  Wiley  and  Sons,  New  York. 


15.  Calculate  an  event  score  for  each  critical  component.  These  scores, 
identified  as  A,y,  are  given  by 

=  A,rA,«  C,- sin  (*„),-  C,  J,  Q,  (11) 

16.  Square  each  component  score.  Accumulate  each  score  and  each 
squared  score  in  bins  reserved  for  these  operations. 

17.  Calculate  the  scores  for  a  total  of  I  similar  sample  events  by 
reiterating  Steps  1-16. 

18.  Applying  the  Law  of  Large  Numbers2,  calculate  an  estimate  of  the 
expected  value  of  kill  for  each  critical  component  as  the  mean  of  the  associ¬ 
ated  scores,  that  is 

=  7  £  x,,  =  x,.  (12A) 

A  bar  is  placed  over  a  quantity  to  identify  its  mean  estimate. 

19.  Calculate  an  estimate  6E{K}}  of  the  standard  deviation  of  the 
estimated  expected  value  of  kill  for  each  component  as10 

<£W'  =  ["A/- 1) -J  ■  (12B) 

Applying  the  Central  Limit  Theorem2,  this  standard  deviation  can  be  inter¬ 
preted  as  implying  that  approximately  68  percent  of  a  large  number  of 
similar  estimates  will  fall  within  one  standard  deviation  of  the  associated 
universe  expected  value. 

Step  19  completes  an  outline  of  an  adjoint  Monte  Carlo  evaluation  of 
the  adjoint  expected-value-of-kill  integral.  An  example  two-dimensional 
problem  is  described  in  the  next  section  of  this  report  and  its  forward-  and 
adjoint-calculated  answers  are  compared. 
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IV  THE  SOLUTION  OF  A  PEDAGOGICAL  EXAMPLE  AIR-DEFENSE- 
END-GAME  PROBLEM  USING  FORWARD  AND  ADJOINT  SAMPLING 
PROCEDURES 


A.  The  Problem-Set  Description 

We  have  devised  a  two-dimensional  problem  (Figure  8)  whose  Monte 
Carlo  solution  illustrates  in  a  simple  manner  the  procedures  outlined  in  the 
preceding  section  of  this  report.  This  problem  is  composed  of  the  following 
parts: 

1.  A  Gun  Coordinate  System  (GCS)  is  defined  for  a  gun  and  radar- 
tracking  system  where  the  muzzle  of  the  gun  is  taken  to  be  the  origin.  The 
location  of  the  gun  muzzle  is  assumed  to  be  precisely  known  with  respect  to 
the  radar  so  that  all  uncertainties  associated  with  locating  the  target  in  the 
GCS  are  due  to  inaccuracies  in  the  radar  fixes.  We  will,  depending  on 
analytical  convenience,  use  either  the  Cartesian  variables  (z.y)  or  the  polar 
variables  (0,s)  (Figure  9)  as  the  coordinates  of  the  bullet  position  vector  rb 
and  the  target-circle-center  position  vector  r,  in  the  GCS. 

2.  A  circular,  two-dimensional  target  of  radius  a  is  moving  in  the  posi¬ 
tive  y  direction  along  the  line  xa  —  2\/36  with  a  constant  velocity  of 
va  =  b/  sec  where  b  is  a  length  (Figure  8).  The  target  center  is  located  at 
( xa  =  2v/36,  yt  =  0)  at  time  t  =  0.  A  Target  Coordinate  System  (TCS)  is 
defined  whose  origin  is  the  the  center  of  the  target  circle  and  whose  vari¬ 
ables  are  taken  to  be  the  polar  coordinates  (<f>,s)  (Figure  10). 

3.  The  radar  detection  and  tracking  of  the  target  is  assumed  to  begin 
at  ya  =  0  and  to  continue  until  a  bullet  is  fired.  These  radar  measured  tar¬ 
get  locations  and  velocities  are  presumed  to  be  used  in  some  unspecified 
manner  to  predict  future  target  locations  which  in  turn  are  used  to  calcu¬ 
late  the  aiming  angle  for  a  time  of  firing  (<,.  In  these  calculations,  the  x- 
coordinate  of  the  target  location  and  the  x-component  of  the  target  velocity 
are  assumed  to  be  precisely  measured  so  that  all  radar  predictions  place  the 
target  center  on  the  line  *,  =  2>/3 b.  Similarly,  during  the  calculation,  the 
acceleration  of  the  target  along  its  path  is  taken  to  be  precisely  zero. 

4.  A  bullet  is  fired  at  some  time  between  /  =  l  sec  and  t  =  2  sec 
where  the  PDF  of  firing  times  is  given  by  the  triangular  distribution  fx(t0) 
(Figure  11a).  The  target  center,  y0>  at  time  of  firing  is  located  at 
(y«)o  =  %  =  v»  to  The  uncertainty  in  the  y-direction  associated  with  the 
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Figure  8.  A  Two-Dimensional  Air-Defense  End-Game  Example  Problem 
(Not  to  Scale) 
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Figure  9.  The  Gun-Coordinate  System  (GCS)  for  the  Example  Problem 
(Not  to  Scale) 
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Figure  11.  The  Probability  Density  Functions  Associated  with  the  Example 
Problem 
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radar  measured  (rm)  target-center  location,  yrm,  is  identified  at  the  time  of 
firing,  /0,  by  the  triangular  conditional  PDF,  /2(yrmiy0)  (Figure  lib),  and  the 
uncertainty  in  the  y-direction  associated  with  the  radar  measured  target 
velocity  vrm  is  described  by  a  similar  triangular  PDF  /3(v,J  (Figure  11c). 

5.  The  aiming  angle  ?0  used  for  a  firing  event  is  derived  from  the  meas¬ 
ured  values  for  target  velocity  and  location  and  the  mean  value,  Fro,  of  the 
bullet  muzzle  velocity. 

6.  The  true  firing  angle  60  due  to  the  various  uncertainties  associated 
with  positioning  and  firing  the  gun  is  predicted  by  the  triangular  condi¬ 
tional  PDF  fM$0)  (Figure  lid). 

7.  The  true  muzzle  velocity  vm  of  the  fired  bullet  has  an  uncertainty 
which  is  predicted  by  the  triangular  PDF  /5(t>m)  (Figure  lie).  The  quantity 
vm  is  the  mean  muzzle  velocity  for  a  large  number  of  similar  firing  events. 

8.  The  bullet  travels  in  a  straight  line  with  a  velocity  which  decreases 
at  a  constant  rate  given  by 

(13A) 

where  tf  is  the  time  that  the  bullet  has  been  in  flight.  The  distance  |rfc|  trav¬ 
eled  by  a  bullet  at  a  time  of  flight  tj  is  derived  from  Equation  12A  as 

N  =  »« </(*  -  o-1  '/*•  (13B) 

9.  Associated  with  the  target  is  a  response  function  Q  which  in  air- 

defense  studies  would  quantify  some  specified  loss  of  capability  by  a  target 
component  due  to  a  hit.  In  this  example  problem,  we  assume  that  the  loss 
of  capability  is  a  function  only  of  the  location  of  a  hit  on  the  target  circle. 
This  response  function  has  the  form,  Q(7[(r4)iy)(r-y,  in  the  GCS  and  the 
simpler  form,  in  the  TCS.  We  identify  the  position  vectors  associ¬ 

ated  with  target  hits  by  affixing  the  subscript  H.  If  QG  and  QT  are  unity  for 
all  hits,  then  the  expected-value  integral  constructed  below  is  merely  the 
expected  value  of  hit  for  the  specified  firing  events. 

Item  9  completes  the  description  of  the  example  problem.  The 
“natural”  formulation  of  the  expected  value  E{Target  Response}  of  target 
response  is  given  by  the  multiple  integral 
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E{Target  Response}  = 


/*  ^#//i(*o)  hiyJiVo)  /a(0  /tf^ol^o)  AK)  X 

<2c((ri)tf>(p«y  ‘fyn  ^0  ^vrm  d Vrm  (14) 

where  rtg  is  the  region  in  the  GCS  space  {toyyrm,vrm,e0,vm)  where  bullets  will 
hit  the  target  circle. 


B.  The  Forward  Monte  Carlo  Estimation  Procedures 


The  multiple  integrals  in  Equation  14  can  easily  be  evaluated  using  a 
“natural'5  Monte  Carlo  procedure  by  simulating  a  large  number  of  firing 
events  using  the  given  PDF’s  and  then  calculating  the  mean  target 
response.  We  outline  below  (Figure  12)  such  a  procedure  for  calculating  the 
expected  value  of  hit  E{H}  by  taking  <?[(rt)H,(r,y  to  be  1  for  all  hits: 

1.  Pick  a  value  for  the  firing  time  of  the  bullet  by  sampling  the  PDF 
/,(*o).  This  sample  value  is  identified  as  (Iq),-  where  » is  the  sample  index  and 
/  is  the  number  of  Monte  Carlo  histories  to  be  used  in  the  calculation.  Cal¬ 
culate  the  true  target  location  at  (^),  as  (%),  =  vt  (*0),. 

2.  Pick  a  value  for  the  radar  measurement  of  the  target  location  at  the 
time  of  firing  by  sampling  the  PDF  /2[yrml(j/0)J-  As  noted  earlier,  this  sample 
value  of  the  measured  target  location  is  always  taken  to  lie  on  the  true 
flight  path  with  an  uncertainty  only  in  the  y-coordinate.  This  sample  value 
is  identified  as  (yrm),. 

3.  Pick  a  value  for  the  measured  target  velocity  by  sampling  the  PDF 
A( *Vm)-  This  sample  value  is  identified  as  (vrm),-  As  noted  earlier,  this  sample 
value  of  the  measured  velocity  is  always  taken  To  be  in  the  positive  y- 
direction. 


4.  Assuming  that  target  locations  at  future  times  are  derived  from  the 
sample  values  (yrm),  and  (rrm)„  calculate  the  aiming  angle  for  which  a  bullet 
having  the  mean  muzzle  velocity  will  hit  the  center  of  the  target.  This 
can  be  accomplished  by  inserting  these  sample  values  into  the  relation 


2 

(1262  -f  (yrw  -f  Vrmtfi* 
Fm  (1-0.1  tj) 


(15) 
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Figure  12.  The  Forward  Monte  Carlo  Estimation  of  the  Expected-Value- 
of-Hit  for  the  Example  Problem 
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derived  from  the  time-of-flight  analytic  right  triangle  (Figure  13),  and  then 
iteratively  locating  the  time  of  flight,  to  a  hit  as  the  fixed  point11  of 
g(tj).  An  aiming  angle  (?0),  is  then  derived  for  the  anticipated  location  of  the 
target  at  the  time  of  hit. 


5.  Pick  a  value  for  the  true  firing  angle  by  sampling  the  conditional 
PDF  /4[Wo),]-  This  sample  value  is  identified  as  (0O),.  As  noted  earlier,  the 
birth  site  of  each  bullet  is  taken  to  be  the  origin  of  the  GCS. 

6.  Pick  a  value  for  the  muzzle  velocity  of  the  fired  bullet  by  sampling 
the  PDF  /5(rm).  This  sample  value  is  identified  as  (t>m),. 


7.  The  determination  of  whether  a  bullet  hits  the  target  is  not  as  easy 
as  might  be  expected  since  both  the  target  circle  and  bullet  are  moving.  In 
the  method  used  here,  we  first  note  that  the  square  of  the  distance  M[t{) 
between  the  bullet  and  the  target-circle  center  is  given  by 

w,)  =  [*,((,)  -  2v^  6]2  +  [y.lf/l  -  Sf.(«))2.  ( 1 6) 

where  the  bullet  time  of  flight  tf  is  related  to  the  target  time  t  by 

t  =  t0+  tj.  (17) 

The  time  of  flight  at  the  distance  of  closest  approach  is  given  by  setting  the 
time  derivative 

d\f[tf) 
dt} 

and  solving  the  resulting  equation  for  its  zero.  Performing  the  indicated 
differentiation  and  grouping  terms  for  convenience,  the  function  g(tf)  is 
derived  such  that 

2v/3  6  -  H[tf)  _ 
tj~  vm  (1  -  0.1  tf)  sin0o  W 

where 


(18A) 


(y»(*/)  -  y«(</)l  (t>m  (i  -  0  2/;)  cosgp  -  t>J 

7  vm  (1  -  0.2fy)  sin0o 


(18C) 


The  function  g(tf)  was  determined  by  trial  to  be  an  acceptable  iteration 
function  for  finding  the  fixed  point  of  Equation  18B  where  the  fixed  point, 
( lj)m ,  is  the  time  of  flight  at  which  the  miss  distance  is  a  minimum. 


Figure  13.  The  Time-of-Flight  Analytic  Right  Triangle  Used  To  Calculate 
Aiming  Angles  (Not  to  Scale) 
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8.  Calculate  the  location  of  the  bullet  and  the  target  at  time  (tf)m. 
Assign  a  score  of  one  to  the  sample  event  if  the  bullet  is  inside  the  target 
circle,  otherwise  assign  a  score  of  zero.  The  score  is  identified  as  Q,. 

9.  Reiterate  Steps  1-8  for  a  total  of  /  histories. 

10.  Calculate  the  expected  value  for  hitting  the  target  as 

£{«}  =  ?*=  4  E  <19A) 

J  —  i 


11.  Calculate  the  standard  deviation  SPH  of  PH  as 


f  E 


\_ 

2 


L  /(/- 1)  J 


(19B) 


C.  The  Adjoint  Monte  Carlo  Estimation  Procedures 


We  will  now  attempt  to  change  some  of  the  variables  of  integration 
used  in  Equation  14  to  obtain  a  form  such  that  bullets  which  hit  the  target 
can  always  be  picked  as  samples. 


To  accomplish  this  objective,  we  first  derive  a  second  analytic  right  tri¬ 
angle  (Figure  14)  which  describes  bullet  intersections  of  the  target  circle. 
We  then  note  that  the  pair  of  independent  equations, 


tan  0Q  = 


2y/3 b  -  a  sin  #  _  2y/3 b  -  a  sin  # 

yt(t)  -  a  cos  4>  vjt  -  a  cos  <f>  ' 


(20A) 


//l -0.1//) 


£ 

{(2v/36  -  a  sin#)2  +  (t>,  t  -  a  cos#)2]2 
//(1-0.K/) 


(20B) 


can  be  derived  which  relate  (0o,vm)  in  the  GCS  to  (#,<)  in  the  TCS.  The 
quantity,  |rt|#,  in  Equation  20B  is  the  distance  traveled  by  a  bullet  to  a  hit. 


We  assume  that  each  realistic  bullet  trajectory  (and  target  maneuver) 
for  a  particular  set  of  values  for  «o,?0)  can  uniquely  determined  by  either 
(0o,vm)  or  (#,<),  so  that  all  points  in  R%  space  (^y^v^O 0,vm)  will  transform 
one-to-one  to  points  in  Rjj  space,  (tQ,yrm,vrn,<l>,l)-  When  this  condition  is 


-46- 


v.y  '..•.'.V.  n .  *  .1  ■  ■  ■  ******* 


v-v 


Figure  14.  The  Time-of-Flight  Analytic  Right  Triangle  Which  Relates 
Target  Hits  in  the  GCS  and  the  TCS  (Not  to  Scale) 


-47’ 


satisfied,  the  expected-value-of-target-response  integral  will  also  be  given  by 
E{  TargtlRtsponse}  = 

/•r't  •f/M  ttPrM  /s(0  AWo)  /s(«m)  M,t)  QAM  d4>  dt  dvrm  dVrn  dt0.  (21) 

We  have  now  described  the  target  response  in  the  target  coordinate  system 
by  the  new  function  Qq{<t>,a).  The  Jacobian  of  the  transformation  is  the 
absolute  value  of  the  determinant  of  the  indicated  partial  derivatives,  that 
is 

d0o  d0o 
d<f>  dt 

dvm  dvm 

d<f>  dt 

The  6q  partial  derivatives  in  the  preceding  determinant  can  be  derived  from 
Equation  20A  to  be 


a  WO  cos  <$>  -f  2\/3  sin  -  a)  sin20o 
d<t>  [ya(/)  -  a  cos  4>)2 

d00  _  va  [2\/3  6  -  a  sin  <f>\  sin20o 
dt  [yjt)  -  a  cos  <t>\2 


(23A) 


(23B) 


Similarly,  the  vm  partial  derivatives  can  be  derived  from  from  Equation  20B 
to  obtain 


and 


dvm  a  (ya(/)  sin  <f>  -  2\/3 b  cos  <t>\ 

W  ~  (  r*)*  t  (1  -0.lt) 


dvm  _  t  v0  [1  -  0.1/]  (ya(t)  -  a  cos  4>]  -  [( r6)^]  [l  -  0.2/] 
81 


I(r6)*]  Z2  (1  -0.1/p 


(23C) 


(23D) 


The  integrals  in  Equation  21  can  also  be  evaluated  by  using  the  Monte 
Carlo  method  where  sample  events  are  now  in  part  determined  by  the  sam¬ 
ple  values  picked  for  the  new  variables  of  integration  <f>  and  t.  Such  a  pro¬ 
cedure,  where  the  first  four  steps  are  identical  to  those  used  earlier  in  the 
“natural”  evaluation,  is  given  below  (Figure  15): 

1.  Pick  a  value  for  the  firing  time  of  the  bullet  by  sampling  the  PDF 
hit o)-  This  sample  value  is  identified  as  (/0),.  Calculate  the  location  (y0),  of 
the  target  center  at  the  sample  firing  time. 
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2.  Pick  a  value  for  the  radar  measurement  of  the  target  location  at  the 
time  of  firing  by  sampling  the  PDF  /i(y„„|(yo).j-  This  sample  value  is 
identified  as  (yrm),  where  i  is  the  sample  index  and  /  is  the  number  of  his¬ 
tories  to  be  used  in  the  calculation. 

3.  Pick  a  value  for  the  measured  target  velocity  by  sampling  the  PDF 
/2(0-  Th  is  sample  value  is  identified  as  (t»rm),-. 

4.  Assuming  that  target  locations  at  future  times  are  derived  from  the 
sample  values  (yrm),  and  (efm)„  calculate  the  aiming  angle  for  which  a  bullet 
having  the  mean  muzzle  velocity  Fm  will  hit  the  center  of  the  target.  This 
can  be  accomplished  by  inserting  these  sample  values  into  the  relation 

[12&2  +  (yrm  +  vmtf)Y  _ 


derived  from  the  time-of-flight  analytic  right  triangle  (Figure  13),  and  then 
iteratively  locating  the_time  of  flight,  ( tf)H ,  to  a  hit  as  the  fixed  point11  of 
g(tf).  An  aiming  angle  (0O),  is  then  derived  for  the  anticipated  location  of  the 
target  at  the  time  of  hit. 

5.  Step  4  is  the  last  step  in  the  adjoint-evaluation  procedure  which  is 
the  same  as  the  forward-evaluation  procedure  described  earlier  in  this  Sec¬ 
tion.  At  this  stage,  the  determination  of  the  i*h  event  score,  X„  could  be 
regarded  as  the  evaluation  of  the  integral, 


\  =  /  fM  I  (*<>),]  4(0  MA  ^  dt , 


(24A) 


where  1S  ^e  reS*on  *n  the  ( <p,t )  subspace  which  includes  all  points 

[Po)»>(yrm)i.(*Vm)><M  contained  in  R%.  If  one  wished,  the  integral  could  be 
evaluated  using  Monte  Carlo  procedures  to  a  high  precision  as 


—  ~ V  ^  U\%)%  I  (^o)J  /sKOn)  Qli4> va)  . 


(24B) 


where  (^B,t„)  are  picked  with  equal  probability  over  Rjt,  and  (0O)„  and  (rm)B 
are  calculated  using  (<*>„,*„).  The  quantities  n  and  N  are  introduced  as  a  new 
running  index  and  the  number  of  associated  samples,  respectively,  and 

A[rtJ\  is  the  area  of  *j),"  However,  applying  the  Law  of  Large  Numbers2, 

an  acceptable  iik  score  can  be  derived  using  a  single  sample  from  the 
integral  in  Equation  24A,  that  is 


■ v>'AV«, 


(24C) 


K  =  [»*;].  AI(#oli  I  («.  l/s«<U]  MM  Qi(M  . 

where  we  have  dropped  the  temporary  index  n  and  reintroduced  the  index 
«.  We  will  describe  the  latter  method  although  the  average  of  five  sample 
values  was  used  in  our  actual  calculations. 

We  would  prefer  to  pick  sample  values  for  with  equal  probability 
at  any  point  in  .,  but  its  boundary  may  present  calculational 
difficulties.  Consequently,  we  will  pick  from  a  larger  region,  probably 
extending  outside  .,  which  presents  no  calculational  difficulty,  and  then 

discard  all  sample  points  ((tol.^Sr which  lie  outside  •  by  the 

simple  expedient  of  setting  to  zero  an  appropriate  factor  which  is  used  later 
to  calculate  the  event  score. 

For  the  sample  values  [( <o),,(yrm)  (OJ.  just  picked  and  the  associated 
calculated  aiming  angle  (0O)„  determine  the  range  (Rt)-  of  values  for  the  tar¬ 
get  time  which  include  all  bullet-target  intercepts  associated  with  non-zero 
values  for  the  PDF's  /4f0oPoM  an<^  fb[(vm)t\-  The  width  of  this  range  is 
identified  as  At,  and  will  be  used  as  a  factor  in  the  calculation  of  the  event 
score. 

6.  Pick  a  sample  value  of  t  with  equal  probability  on  the  sampling 
range  defined  in  Step  5.  This  sample  time  is  identified  as 

7.  Locate  the  target  circle  at  f,  and  calculate  the  range,  (/?^)„  (Figure 
16)  subtended  by  that  part  of  the  target  circle  which  can  be  impacted  by  a 
bullet.  The  width  of  this  range  is  identified  as  A<£,  and  will  be  used  as  a  fac¬ 
tor  in  the  calculation  of  the  event  score. 

8.  Pick  a  value  for  <f>  over  the  range  calculated  in  Step  7.  This  sample 
value  is  identified  as  <£,. 

9.  Calculate  in  the  GCS  the  location  of  the  target  center  at  time 
Calculate  the  location  of  the  point  on  the  target  circle  associated  with  (<£,-, f,). 
The  Cartesian  coordinates  of  this  point  are  identified  as  [{%),,( yw),j.  Since  all 
bullet  trajectories  intersect  the  target  circle,  the  quantity,  Qrl(4>,),a], 
identified  later  as  is  set  to  1  for  all  events. 

10.  Calculate  the  value  of  90  associated  with  [(*#)„( yw),].  This  value  is 
identified  as  (0O),. 
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11.  Calculate  the  value  of  the  PDF  /4((0o).-l(^o).,l-  This  value  is  identified 
as  p,{0o)  and  will  be  used  as  a  factor  in  the  calculation  of  the  event  score.  A 
zero  value  implies  that  the  sample  firing  event  lies  outside  /?£  and  should  be 
discarded.  Jump  to  Step  17  if  p,{60)  =  0;  otherwise,  proceed  to  Step  12. 

12.  Calculate  the  value  of  vm  such  that  the  sample  firing  will  deliver  a 
bullet  at  [(!//),,( y//),]  at  time  f,.  This  value  is  identified  as  (vm),-. 

13.  Calculate  the  value  of  the  PDF  /5((vm).i  This  value  is  identified  as 
pt(vj  and  will  be  used  as  a  factor  in  the  calculation  of  the  event  score.  A 
zero  value  implies  that  the  sample  firing  event  lies  outside  R]j  and  should  be 
discarded.  Proceed  to  Step  17  if  p,{vm)  =  0;  otherwise,  proceed  to  Step  14. 

14.  Calculate  the  Jacobian  for  the  history.  This  quantity  is  identified 
as  J{. 

15.  Calculate  the  event  score  X,  as 

X,  =  At,  A0,  p,{0o)  p,{vm)  J,  Qi.  (25) 

Accumulate  the  score  in  a  bin  reserved  for  this  operation. 

16.  Square  the  event  score  X,  and  accumulate  the  result  in  another  bin 
reserved  for  this  operation. 

17.  Reiterate  Steps  1-16  for  a  total  of  /  histories. 

18.  Applying  the  Law  of  Large  Numbers2,  calculate  an  estimate  of  the 
probability  of  hitting  the  target  as 

Ph  =  4  E  (26A) 

1  1-1 


19.  Calculate  an  estimate  6PH  of  the  standard  deviation  of  the  hit- 
probability  estimate  as10 


E 

-  1 _ 

1(1-  1) 


2 


(26B) 


Step  19  completes  a  general  outline  of  the  Monte  Carlo  evaluation  of 
the  multiple  integrals  used  to  formulate  the  expected  value  of  hit  for  the 
example  problem.  We  will  discuss  the  answers  to  this  problem  below  in 
Section  IV. B. 
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C.  The  Problem  Results 


The  FORTRAN  computer  programs  FORWARD  and  BACK  were  con¬ 
structed  and  used  to  calculate  in  turn  forward  and  adjoint  estimates  of  the 
expected-value-of-hit  integrals  given  earlier  in  this  section.  A  sufficient 
number  of  sample  histories  (5000)  to  obtain  an  acceptable  standard  devia¬ 
tion  for  the  set  of  problems  described  below  was  used.  In  fact,  the  stan¬ 
dard  deviations  obtained  using  this  large  number  of  sample  histories  are 
generally  much  smaller  than  those  obtained  in  most  air-defense  end-game 
calculations.  The  calculational  time  needed  by  BACK  was  estimated  to  be 
approximately  twice  that  of  FORWARD  for  an  equal  number  of  sample 
histories. 


Using  each  computer  program,  we  calculated  the  expected  value  of  hit 
for  a  set  of  targets  whose  radii  varied  systematically  from  a  ==  0.01056  to 
fl  =  0.00016  where  6  is  taken  to  be  1000  (arbitrary  units).  The  target  velo¬ 
city,  v#,  was  taken  to  be  b/sec,  and  the  mean  bullet  velocity  vb,  was  taken 


to  be  3b/sec.  The  fractional  standard  deviation  percentage 


6P  u 

100  -=l 
Ph 


L  J 

remained  approximately  constant  for  the  set  of  answers  calculated  using 
BACK;  while  the  fractional  standard  deviation  percentage  calculated  using 
FORWARD  increased  with  decreasing  target-circle  radius  (and  conse¬ 
quently  decreasing  PH).  This  result  was  expected  since  a  large  fraction  of 
the  scores  in  the  adjoint  calculation  have  intermediate  values  while  only  O’s 
and  l’s  are  scored  in  the  forward  mode.  These  answers  are  tabulated  below 
in  Table  1. 


The  desirability  of  using  the  adjoint  method  to  calculate  the  expected 
value  of  hit  where  this  value  is  small  can  be  demonstrated  by  graphing  the 

bPH 

weighted  standard  deviations  W-=^~  vs.  a  for  the  twx>  calculational  methods 

P H 

(Figure  17).  The  weight  (calculational-time  ratio)  is  taken  to  be  2  for  the 
BACK  computer  program  and  1  for  the  FORWARD  computer  program. 
The  crossover  point  at  a  =  3.4  defines  the  point  at  which  one  should  con¬ 
sider  using  adjoint  sampling  procedures. 


ft 


ADJOINT  CALCULATION  TABLE  A  .THE  FORWARD  CALCULATION 


0 

Ph 

SPH 

IOOSPh/Ph 

10.5 

0,00688 

1,116 

8.5 

0.50660 

0.00707 

1,396 

6.5 

0,40760 

0.00695 

1.705 

4,5 

0.27840  • 

0,00634 

2.277 

2,5 

0.16540 

0,00525 

3.177 

2.0 

0.12560 

0,00469 

3,732 

1.5 

0.10420 

0,00432 

4.147 

1  .0 

0.06420 

0.00347 

5.400 

0.5 

0,03400 

0.00256 

7,539 

0.4 

0.02520 

0,00222 

8.797 

0.3 

0,02140 

0.00205 

9,564 

0.2 

0,01320 

0.00161 

12.229 

_  0.1 

0.00560 

0,00106 

18,847  _ 

~~  10.5 

0,62371 

0.00926 

1,484 

8.5 

0,53351 

0.00753 

1,412 

6.5 

0.41915 

0,00584 

1,392 

4.5 

0,29440 

0.00400 

1,359 

2.5 

0.16662 

0,00223 

1,339 

2.0 

0,13252 

0.00180 

1.356 

2.5 

0.09870 

0.00135 

1,372 

1.0 

0.06539 

0,00090 

1,372 

0.5 

0,03222 

0,00044 

1.378 

0.4 

0,02532 

0.00035 

1.396 

0.3 

0,01977 

0,00027 

1.377 

0.2 

0,01315 

0,00018 

1,351 

0.1 

0,00640 

0,00009 

1,382 

Table  1.  The  Example-Problem  Answers  Calculated  Using  Both  Forward 
and  Adjoint  Procedures 
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Comparison  of  the  Calculational  Efficiencies  of  the  Forward 
t  Monte  Carlo  Procedures 


V.  CONCLUSIONS  AND  PROGNOSIS 


We  have  constructed  an  integral  representation  of  the  expected  value  of 
kill  of  a  critical  component  in  an  aircraft  by  a  bullet  fired  from  a  gun 
located  on  the  ground,  and  have  outlined  a  Monte  Carlo  evaluation  of  this 
integral  where  sample  events  are  picked  from  their  “natural”  distributions. 
We  then  changed  the  variables  of  integration  (0o,iJ>Q,c4)  (firing  angles  and 
air-drag  coefficient)  used  in  the  “natural”  form  of  the  expected-value-of-kill 
integral  to  (polar  and  azimuth  angle  in  the  target  spherical  coordinate 
system  plus  the  target  time),  and  outlined  the  Monte  Carlo  evaluation  of 
the  transformed  integral.  The  necessary  probability  density  functions  for 
picking  sample  values  for  the  variables  of  integration  are  identified  and  the 
bullet  flyout  is  described  by  equations  given  in  Reference  1.  However,  the 
outlined  procedures  are  not  restricted  to  any  particular  description  of  bullet 
flyout,  but  should  be  applicable  to  other  methods  of  predicting  the  bullet 
flyout. 

A  quantitative  comparison  of  the  calculational  efficiency  of  the  two 
methods  for  a  particular  air-defense  problem  can  only  be  obtained  by  trial 
calculations.  In  conducting  a  calculation,  the  effort  of  promenading  a  bullet 
near  the  target  to  obtain  the  interaction  point  would  be  expected  to  com¬ 
pare  with  the  effort  needed  to  calculate  the  values  of  the  changed  variables 
of  integration  by  solving  three  simultaneous  equations.  However,  as  illus¬ 
trated  for  the  set  of  example  problems  in  Section  IV,  the  fractional  stan¬ 
dard  deviation  of  a  hit  (or  kill)  probability  estimate  would  be  expected  to 
increase  with  decreasing  hit  probability  when  forward  sampling  is  used,  and 
to  remain  approximately  constant  with  decreasing  hit  probability  when 
adjoint  sampling  is  used.  If  the  forward  method  is  more  efficient  for  prob¬ 
lems  having  a  hit  probability  near  one,  then  a  crossover  exists  at  some  hit 
probability  beyond  which  the  adjoint  method  would  become  the  more 
efficient  method.  To  determine  when  adjoint  sampling  is  the  more  efficient 
method,  the  answers  should  be  calculated  for  a  set  of  representative  air- 
defense-end-game  problems  similar  to  the  example  set  described  in  Section 
IV. 


Our  complete  omission  of  geometrical  detail  in  the  target  when  describ¬ 
ing  the  interaction  of  bullets  or  shrapnel  with  its  critical  components  has 
been  deliberate.  A  combinatorial  geometry1*  description  of  military 


**W.  Guber,  R.  Nagel,  R.  Goldstein,  P.S.  Mittelman,  and  M.  Kalos,  “A  Geometric 
Description  Technique  for  Computer  Analysis  of  Both  the  Nuclear  and  Conventional 
Vulnerability  of  Armored  Military  Vehicles,”  Mathematical  Applications  Group,  Inc., 
Report  No.  MAGI-6701,  August  1967. 
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vehicles  can  be  constructed  to  a  high  precision,18  14  but  the  validity  of  the 
present  procedures16  being  used  to  transport  penetrators  into  surface  tar¬ 
gets  and  the  prediction6  of  the  incapacitation  of  an  impact  with  a  critical 
component  have  not,  to  the  knowledge  of  the  author,  been  rigorously 
verified.  In  the  meantime,  we  expect  that  greatly  simplified  target  descrip¬ 
tions,  such  as  those  used  by  the  current  MGEM1  computer  program  pro¬ 
vide  all  the  target  detail  that  is  warranted. 

A  causal  description  of  the  detonation  of  a  PX  bullet  and  the  transport 
of  the  resulting  shrapnel  into  the  target  aircraft  to  possible  impacts  with 
critical  components  would  require  a  prohibitive  effort.  A  stochastic 
approach  has  been  suggested  in  Reference  8  whereby  the  behavior  of  a  large 
aggregate  of  similar  ballistic  events  may  be  adequately  describable  by  pro¬ 
bability  density  functions  of  properly  selected  random  variables.  The 
existence  of  such  a  set  of  random  variables  has  not  been  verified,  but  the 
adequacy  criterion  would  not  be  expected  to  be  very  severe  if  the  majority 
of  the  metal  fragments  impacting  a  critical  component  have  perforated  only 
one  or  two  solid  barriers.  In  any  event,  the  adequate  description  of  such 
ballistic  effects  is  essential  to  vulnerability  analysises,  and  a  serious  attempt 
should  be  made  to  develop  and  verify  such  a  predictive  model. 

If  such  a  stochastic  description  of  ballistic  events  is  developed  and 
verified,  then  calculational  methods  developed  for  radiation  transport1® 17 


13L.  Bain,  and  M.  Reisinger,  “The  GIFT  Code  User  Manual;  Vol.  I,  Introduction  and 
Input  Requirements,”  Ballistic  Research  Laboratory  Technical  Report  No.  1902,  July 
1975,  (UNCLASSIFIED). 

l4G.  Kuehl,  M.  Reisinger,  and  L.  Bain,  “The  GIFT  Code  User  Manual;  Vol.  II,  The 
Output  Options,”  Ballistic  Research  Laboratory  Technical  Report  No.  02189,  September 
1979,  (UNCLASSIFIED). 

,6T.S.  Hafer  and  AS.  Hafer,  “Vulnerability  Analysis  for  Surface  Targets  (VAST):  An 
Internal  Point  Burst  Model,”  Ballistic  Research  Laboratory  Technical  Report  No.  02154, 
April  1979,  (UNCLASSIFIED). 

1SJ.  Spanier  and  E.M.  Gelbard,  Monte  Carlo  Principles  and  Neutron  Transport 
Problems.  Addison-Wesley  Publishing  Company,  Reading,  Massachusetts,  Menlo  Park, 
California,  London,  Don  Mills,  Ontario. 

l7W.C.  Roesch,  “Mathematical  Theory  of  Radiation  Fields,”  Basic  Concepts  of 
Dosimetry.  F.H.  Attix  and  W.C.  Roesch,  Editors,  Volume  I,  Fundamentals,  Second 
Edition,  Academic  Press,  New  York,  San  Francisco,  and  London,  1968. 
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may  be  advantageously  applied  to  ballistic  vulnerability  studies.  For  exam¬ 
ple,  looking  at  the  shrapnel  generated  during  the  detonation  of  a  PX  bul!  t, 
we  would  expect  that  the  formulation  of  the  expected-value-of-kill  integral 
would  permit  the  calculation  of  the  adjoint  fragment  fluence,  $r0,w0,b0), 
over  the  fuzing  region,  R},  for  each  critical  component  in  the  aircraft.  The 
kill  probability  per  fragment,  PK,  would  then  be  given  as 

PK  =  Ilf 5/r0,w0,b0)  4>(r0,w0,b0)  A>0rfw0<ft>0,  (27 A) 

*/ 

Where  Sf(r0,w0,b0)  is  the  probability  density  of  shrapnel  fragments  produced 
during  a  fuzing  action.  The  quantities,  (r0,w0),  are  the  position  vector  and 
direction  of  motion,  respectively,  of  a  fragment  in  a  target  coordinate  sys¬ 
tem,  and  b0  is  a  generalized  set  of  variables  which  can  be  used  to  character¬ 
ize  fragments.  In  practice,  0(ro,wo,bo)  would  be  evaluated  in  bins  over  the 
fuzing  region,  Rf,  during  a  one-time  calculation,  and  PK  would  then  be  cal¬ 
culated  as  the  sum  of  the  products,  (S^r), $■.  that  is, 

fy=  E  (%*.••  (27B) 

I  -  1 

The  quantity  i  is  the  running  index  over  bins  in  the  fuzing  region  about  the 
target,  and  /  is  the  total  number  of  such  bins. 

The  advantage  of  this  method  is  that  the  adjoint  fluence  can  be 
evaluated  over  the  fuzing  region  of  a  plane  in  a  one-time  calculation  and 
then  be  reused  for  any  number  of  specific  problems.  A  possible  disadvan¬ 
tage  in  using  the  adjoint  fluence  is  that  its  evaluation  may  be  laborious  and 
the  dimensionality  of  the  mesh  needed  for  its  description  over  the  fuzing 
region  may  be  large.  In  any  event,  the  feasibility  of  such  a  procedure 
should  be  investigated. 
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APPENDIX 


THE  MONTE  CARLO  EVALUATION  OF  EXPECTED- VALUE  INTEGRALS 
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THE  MONTE  CARLO  EVALUATION  OF  EXPECTED-VALUE  INTEGRALS 

The  general  features  of  a  Monte  Carlo  evaluation  of  the  expected- 
value-of-kill  integrals  used  to  describe  the  efficiency  of  air-defense  gun  sys¬ 
tems  are  illustrated  below  by  the  triple  integral 

E{  Q(X,  Y,Z)}=f  /  /  $z,y,z)  Q{it,v,w)  dz  iy  dz  = 

Rrfi 

/  /  /  S{x,y,z)  Q[v(z,y,z),^z,y,z),n{z,y,z)\  dz  dy  dz.  (Al) 

The  quantity  5(z,y,2)  is  the  joint  probability  density  of  values  for  the  set  of 
random  variables  (X,Y,Z)  in  a  region  of  interest  RSfl  of  zyz  space,  and  Q(u,v,w) 
is  any  well-behaved  function  with  a  finite  maximum  value.  In  certain  physi¬ 
cal  problems,  S(z,y,z)  is  often  identified  as  the  source  term,  and  Q{u,v,w )  is 
identified  as  a  pay-off  function  or  detector-response  function.  For  example, 
S(jc,y,2)  could  describe  the  probability  density  of  particles  generated  by  some 
physical  process,  and  Q(u,v,w)  could  describe  the  interaction  of  these  parti¬ 
cles  with  a  detector  after  their  passage  through  an  interactive  medium.  In 
our  air-defense  end  games,  the  source  term  is  taken  to  include  the  uncer¬ 
tainties  associated  with  detecting  and  tracking  a  target,  aiming  the  gun 
and  firing  the  bullet,  and  transporting  the  bullet  toward  the  target;  and  the 
pay-off  function  is  the  kill  probability  of  a  critical  component  in  the  target 
by  the  bullet.  The  source  term  will  have  its  largest  values  in  some  known 
region  and  will  become  vanishingly  small  at  known  distances  from  this 
region.  Therefore,  the  region  RS9,  is  generally  taken  to  include  all  points  in 
zyz  space  whose  inclusion  in  an  expected-value  integral  will  make  a 
significant  contribution  to  the  expected  value.  A  source  term  is  generally 
assumed  to  be  normalized  to  unity  unless  otherwise  noted,  that  is 

/  /  /  S(z,y,z)  dz  dy  dz  =  1.  (A2) 

Rzyi 

According  to  Kaplan16  the  functions 

u  =  u(z,y,2),  (A3) 

v  =  t {z,y,z),  (A4) 

and 


l6W.  Kaplan,  Advanced  Calculus.  Addison-Wesley  Publishing  Company,  Inc., 
Cambridge  42,  Mass,  1953. 


w  =  ti(z,y,z),  (A5) 

are  assumed  to  be  defined  and  continuous  in  the  region  RtJt  of  xyz  space. 

The  corresponding  points  lie  in  the  region  Rutm  of  ww  space,  and  it  is 

assumed  that  the  inverse  relations 

x=  x(u,v,w),  (A6) 

y=y(u,v,w),  (A7) 

and 

z=z(u,v,w),  (A8) 

are  defined  and  continuous  in  the  region  R„m  of  uvw  space,  so  that  the 

correspondence  between  and  Rnt  is  one-to-one.  As  used  in  Equation 
Al,  Equations  A3-A5  causally  describe  some  physical  process.  In  our  parti¬ 
cle  example,  (X,Y,Z)  are  taken  to  quantify  the  essential  features  of  particles 
at  their  generation  (birth  particles),  and  ( U,V,W)  describe  these  same  parti¬ 
cles  (residual  particles)  at  their  interaction  with  the  detector  after  passing 
through  an  interactive  medium. 

In  practice,  the  expected-value  integrals  used  in  air-defense  studies  are 
often  of  much  higher  dimensionality  than  that  used  in  the  example  triple 
integral.  In  such  cases  where  a  confidence  level  of  b%  or  so  is  acceptable, 
the  Monte  Carlo  method  is  quite  often  the  most  efficient  method  of  evalua¬ 
tion  since  convergence  toward  the  true  value  always  goes  as  l/v^V  where  N 
is  the  number  of  sample  histories.  An  outline  of  a  general  Monte  Carlo 
evaluation  of  the  integral  used  in  Equation  1  is  given  below: 

1.  Pick  a  set  of  values  for  (x,y,z)  by  sampling  S(x,y,z).  This  set  of  values 

is  identified  as  where  t  is  the  running  index  over  sample  events. 

2.  Calculate  the  associated  values  of  («,»,»)  by  using  Equations  A3-A5. 

This  set  of  values  is  identified  as  In  a  physical  problem,  this 

operation  is  often  interpreted  as  the  simulation  of  a  physical  event  such  as 
the  transport  of  a  particle  through  an  interactive  medium. 

3.  Calculate  the  score  for  the  event  as  This  quantity  is 

identified  as  <?,.  - — 

4.  Accumulate  Q{  in  a  bin  reserved  for  this  operation.  Square  the  score 
and  similarly  accumulate  the  squared  score  in  a  second  bin. 
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5.  Reiterate  Steps  1-4  until  a  total  of  /sample  events  have  been  simu¬ 
lated. 


6.  Applying  the  Law  of  Large  Numbers*,  calculate  an  estimate  of  the 
triple  integral  in  Equation  A1  as  the  mean  of  the  sample  scores,  that  is 

b(<XU,VM=  Q  =  4se.'  (A9) 

1  i—  1 

The  mean  estimate  of  an  integral  is  identified  by  placing  a  bar  over  its 
symbol. 


7.  Calculate  an  estimate  of  the  confidence  level  of  the  preceding  esti¬ 
mate  as  the  standard  deviation  of  the  mean  estimate,  that  is 


6E{Q(UXW)}  =  &Q  = 


i—i 


4/-1) 


(A10) 


Applying  the  Central  Limit  Theorem*,  this  estimate  of  the  confidence  level 
of  the  calculation  can  be  interpreted  as  implying  that  approximately  68% 
of  a  large  number  of  similar  estimates  will  fall  within  one  standard  devia¬ 
tion  of  the  true  value  of  the  triple  integral. 


Step  7  completes  the  outline  of  a  general  solution  of  the  triple  integral 
in  Equation  Al.  This  triple  integral  can  be  reformulated  as 


E{Q(VXW)}  =  CE{S(XXZ)}  = 

m 

c  J  /  dv  dv  dw  (All) 

D  t 


by  changing  the  variables  of  integration  from  ( x,y,z )  to  [v,v,w).  The  quantity 
J[u,v,w)  is  the  Jacobian  of  the  transformation,  that  is 


Av,v,w) 


dz,dy,dz 

du,dv,dw 


(A12) 


The  quantity  <7  is  used  to  normalize  the  function  Q(u,v, u>)  to  unity  over  this 
region,  that  is 


/  /  /  Qiu>v>w)  dv  dv  dw  =  C.  (A13) 
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The  transformed  triple  integral  of  Equation  All  can  also  be  evaluated 
using  the  Monte  Carlo  method  by  interchanging  the  roles  played  by  the 

functions  S(x,y,z )  and  in  the  Monte  Carlo  simulation  procedure. 

G 

Such  a  procedure  is  composed  of  the  following  steps: 

1.  Pick  a  set  of  values  for  («,«,«>)  by  sampling  over  the  region 

This  set  of  values  is  identified  as  where  i  is  the  running  index 

over  sample  events. 

2.  Calculate  the  associated  values  of  {z,p,z)  by  using  the  set  of  inverse 

equations  derived  from  Equations  A3-A5.  This  set  of  values  is  identified  as 
(*;>»;>*.)•  a  physical  problem,  this  operation  can  often  be  interpreted  as 
the  adjoint  simulation  of  a  physical  event  such  as  the  transport  of  a  parti¬ 
cle  from  a  detector  back  to  its  origin.  — 

3.  Calculate  the  score  for  the  event  as  CJ[ ti„ v#) u>,) S( x„ *,) .  This  quantity 
is  identified  as  CJ.S,  or  CX,. 

4.  Accumulate  CX,  in  a  bin  reserved  for  this  operation.  Square  the 
score  and  similarly  accumulate  the  squared  score  in  a  second  bin. 

5.  Reiterate  Steps  1-4  until  a  total  of  /  sample  events  have  been  simu¬ 
lated. 

6.  Applying  the  Law  of  Large  Numbers*,  calculate  an  estimate  of  the 
triple  integral  in  Equation  All  (and  Equation  Al)  as  the  mean  of  the  sam¬ 
ple  scores,  that  is 

5MW  =  T  E  4  s.  =  c*  (AH) 

1  i—  1 

where 

x  =  E  J,S,  (A15) 

I  —  1 

7.  Calculate  an  estimate  of  the  confidence  level  of  the  preceding  esti¬ 
mate  as  the  standard  deviation  of  the  mean  estimate,  that  is 

£  C2X,*-/C2X*  1 

eE{<XV,V,W))  =  iii - -  .  (A16) 
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Applying  the  Central  Limit  Theorem',  this  sstim.tr  of  the  confidence  lerel 
of  the  calculation  can  be  interpreted  an  implying  that  approximately  68 /0 
of  a  large  number  of  similar  estimates  will  fall  within  one  standard  devia- 
tion  of  the  true  value  of  the  triple  integral. 


i 


DISTRIBUTION  LIST 


No.  of 

Copies  Organization 


No.  of 

Copies  Organization 


2  Administrator 

Defense  Technical  Info  Center 
ATTN:  DTIC-DDA 
Cameron  Station 
Alexandria,  VA  22304-6lL5 

1  HQPA 

DAMA-ART-M 

Washington,  DC  20310 

1  Commander 

Armament  R&P  Center 
US  Army  MICOM 
ATTN:  SMCAR-TSS 
Dover,  NJ  07801 

1  Commander 

Armament  R&D  Center 
US  Army  MICOM 
ATTN :  SMCAR-TDC 
Dover,  NJ  07801 

1  Director 

Benet  Weapons  Laboratory 
Armament  R&D  Center 
US  Army  AMCCOM 
ATTN:  SMCAR-LCB-TL 
Watervliet,  NY  12189 

1  Commander 

US  Army  Armament,  Munitions  & 
Chemical  Command 
ATTN :  SMCAR-ESP-L 
Rock  Island,  IL  61299 

1  Commander 

US  Army  Aviation  Research  & 
Development  Laboratory 
Ames  Research  Center 
Moffett  Field,  CA  9^035 

1  Commander 

US  Army  Communications  - 
Electronics  Command 
ATTN :  AMSEL-ED 
Fort  Monmouth,  NJ  07703 


1  Commander 

ERADC0M  Technical  Library 
ATTN:  DELSD-L  (Reports  Section) 

Fort  Monmouth,  NJ  07703-5301 

1  Commander 

US  Army  Missile  Command  Research 
Division  &  Engineering  Center 
ATTN:  AMSMI-RD 
Redstone  Arsenal,  AL  35898 

1  Commander 

US  Army  Missile  &  Space  Intelli¬ 
gence  Center 
ATTN:  AMSMI-YDL 
Redstone  Arsenal,  AL  35898 

1  Commander 

US  Army  Tank  Automotive  Command 
ATTN:  AMSTA-TSL 
Warren,  MI  48397-5000 

1  Director 

US  Army  TRADOC  Systems  Analysis 
Activity 
ATTN:  ATAA-SL 

White  Sands  Missile  Range,  NM  88002 

1  Commandant 

US  Army  Infantry  School 
ATTN:  ATSH-CD-CS0-0R 
Fort  Benning,  GA  31905 

1  Commander 

US  Army  Development  &  Employment 
Agency 

ATTN :  M0DE-TED-SAB 
Fort  Levis,  WA  98^33 


1  Air  Force  Armament  Laboratory 
ATTN :  AFATL/DLODL 
Eglin  AFB,  FL  325^2-5000 


DISTRIBUTION  LIST 


No.  of 

Copies  Organization 
1  Director 

U.S.  Army  Air  Mobility  Research 

and  Development  Laboratory  < 

Ames  Research  Center  I 

Moffett  Field,  CA  94035  j 

10  Central  Intelligence  Agency 

Office  of  Central  Reference 
Dissemination  Branch 

Room  GE-47  HQS  1 

Washington,  D.C.  20502 

Aberdeen  Proving  Ground,  MD 

I 

Dir,  USAMSAA 
ATTN:  AMXSY-D 

AMXSY-MP,  H.  Cohen  I 

Cdr ,  USATECOM 
ATTN:  AMSTE-TO-F 

Cdr,  CRDEC,  AMCCOM 

ATTN:  SMCCR-RSP-A  i 

SMCCR-MU 
SMCCR-SPS-IL 


72 


USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the 
reports  it  publishes.  Your  comments/answers  to  the  items/questions  below  will 
aid  us  in  our  efforts. 

1.  BRL  Report  Number _ Date  of  Report 

2.  Date  Report  Received _ 


3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or 
other  area  of  interest  for  which  the  report  will  be  used.) 


4.  How  specifically,  is  the  report  being  used?  (Information  source,  design 
data,  procedure,  source  of  ideas,  etc.) _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far 
as  man-hours  or  dollars  saved,  operating  costs  avoided  or  efficiencies  achieved, 
etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future 
reports?  (Indicate  changes  to  organization,  technical  content,  format,  etc.) 


CURRENT 

ADDRESS 


Organization 


Address 


City,  State,  Zip 


7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the 
New  or  Correct  Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 


OLD 

ADDRESS 


Organization 


Address 


City,  State,  Zip 


(Remove  this  sheet  along  the  perforation,  fold  as  indicated,  staple  or  tape 
closed,  and  mail . ) 


w; 


w 


